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ABSTRACT 


Expendable  bathymetric  temperature  ( XBT)  data  taken  from 
an  anticyclonic  meander  crest  within  the  Gulf  Stream  (Hummon 
et  al  1991)  is  analysed  by  looking  at  the  empirical  vertical 
structure.  The  ensemble  averaged  data  is  formed  into  a 
projection  matrix  that  compares  the  value  of  the  temperature 
at  one  depth  with  the  temperature  at  a  second  depth.  The  data 
is  smoothed  with  the  correlation  analysis  being  performed  at 
10  metre  intervals  from  5  metres  to  a  depth  of  800  metres. 
The  first  four,  or  principle,  EOFs  of  the  projection  matrix 
are  computed  and  the  modal  amplitudes  for  each  XBT  determined. 
Using  objective  analysis  the  modal  amplitudes  are  interpolated 
onto  a  specified  grid.  Synthetic  XBTs  are  then  reconstructed 
at  the  grid  positions  using  the  interpolated  grid  modal 
amplitude  values.  A  measure  of  the  error  variance  at  each  grid 
point  is  determined.  The  objective  analysis  is  repeated  using 
successively  fewer  XBTs  from  the  data  set,  until  the  resulting 
error  in  the  interpolated  XBTs  at  the  grid  points  become 


unacceptable. 
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I  INTRODUCTION 


From  a  military  standpoint,  to  carry  out  a  successful 
range  prediction  against  a  surface  or  subsurface  unit,  for 
either  passive  or  active  SONAR,  it  is  necessary  to  have  access 
to  the  most  recent  vertical  temperature  profile  that  is 
available  for  the  area  of  interest.  If  several  XBTs  are  taken 
at  different  positions  and  at  different  times  within  a 
region,  what  is  the  optimal  vertical  profile  at  some  arbitrary 
point  of  significance  within  the  region  based  upon  this 
collected  data?  Or,  perhaps  if  multiple  units  are  on  task, 
each  taking  their  own  XBTs,  what  is  the  optimal  interpolation 
of  the  water  condition  at  some  point  between  the  units?  The 
development  of  range  dependent  prediction  models  makes  a 
knowledge  of  the  water  conditions  between  a  unit  and  its 
target  even  more  crucial.  Thus  the  ability  to  be  able  to 
empirically  assess  the  vertical  water  conditions  at  any  point 
within  a  target  region,  and  to  obtain  valid  and  useful 
information,  is  of  considerable  importance. 

From  a  purely  scientific  basis,  it  would  be  of  benefit  to 
know  the  approximate  number  of  vertical  profiles  that  need  to 
be  obtained  before  a  comprehensive  analysis  of  a  given  area 
could  be  achieved.  Similarly,  some  measure  of  the  optimal 
spacing  between  XBTs  would  be  of  value  for  planning  and  the 
economic  use  of  valuable  assets  and  time. 
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Recent  developments  in  satellite  technology  now  allow  the 
determination  of  the  subsurface  vertical  structure  by 
measuring  the  dynamic  height  of  the  ocean  using  altimetry 
(Carnes  et  al  1990) .  But,  how  many  readings  need  to  be  taken 
for  a  given  region  of  the  ocean  before  a  realistic 
interpretation  can  be  constructed?  Additionally,  if  gaps 
exist  within  the  data  collection,  how  much  error  will  exist  in 
interpolating  data  void  areas?  If  sufficient  remote  readings 
could  randomly  be  taken  in  and  around  a  given  feature,  how 
many  readings  would  be  required  before  the  feature's  vertical 
temperature  structure  can  be  adequately  reproduced? 

The  ultimate  goal  of  this  study  is  to  find  out  how  few 
XBTs  are  required  before  an  adequate  vertical  temperature 
profile  can  be  compiled  within  a  Gulf  Stream  meander. 

The  feature  analyzed  in  this  study  is  the  warm  side  of  a 
Gulf  Stream  meander  that  was  identified  and  rigorously  sampled 
during  1988  (Hummon  et  ai  91) .  It  is  anticipated  that  a 
study  of  this  type  conducted  in  this  particular  area  will  be 
of  general  use,  and  give  an  indication  of  the  number  of  XBTs 
that  need  to  be  deployed  before  an  adequate  interpolation  can 
be  made  as  to  the  underlying  water  structure. 

Carter  and  Robinson  (1987)  considered  empirically  the 
effects  of  reducing  the  size  of  an  original  data  set  upon 
the  value  of  the  contour  maps  that  were  produced.  They 
considered  the  depth  of  the  15  degree  Celsius  isotherm.  The 
data  were  taken  during  the  POLYMODE  experiment  and  consisted 
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of  443  XBTs .  The  depth  of  the  15  degree  isotherm  was  extracted 
from  each  of  these  XBTs  and  optimally  interpolated  onto  a 
regular  grid.  The  results  of  the  interpolation  are  shown  at 
the  top  left  of  Figure  1  with  the  associated  amount  of  error, 
for  any  location,  depicted  in  the  top  right  of  Figure  1.  The 
data  set  was  then  halved  and  the  analysis  repeated  (the  middle 
two  pictures)  ,  with  the  effect  that  the  new  results  showed 
very  little  change  in  the  error  field.  However,  as  the  last 
two  diagrams  show,  by  the  time  only  a  quarter  of  the  data  is 
included  the  error  in  the  analysis  field  has  grown 
considerably,  to  the  point  where  the  analysis  is  unacceptable 
for  practical  purposes.  The  study  indicates  that  in  order 
to  survey  the  given  area  of  the  ocean  it  would  have  been 
sufficient  to  have  launched  half  of  the  XRTs  that  were 
actually  launched  without  any  serious  decrease  in  the  quality 
of  the  15  degree  isotherm  map  that  was  produced.  Additionally, 
the  interpolation  procedure  they  employed  gives  an  explicit 
statement  of  the  amount  of  error  involved  in  the 

reconstruction  at  any  location  within  the  region  thereby, 
giving  an  unequivocal  statement  about  its  usefulness,  or 
otherwise,  to  a  future  user. 

The  object  of  this  study  is  to  consider  the  reduction  of 
data  problem  in  an  objective  analysis  procedure  more 
rigorously.  By  reducing  a  data  set  repeatedly  by  one 
observation  until  the  resulting  error  in  the  reconstructed 
vertical  temperature  profile  becomes  unacceptable,  the  minimum 
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MIN*538  MAX* 707  MIN'O.M  MAX*1.00 

C 


Figure  1  The  result  ot  an  objective  analysis  of  the 
depth  of  the  15  C  isotherm  using  different  amounts  of 
data,  for  a  six  degree  square  centred  at  70  W  29  N.  25  m 
contour  interval  for  analysis,  o.25  contour  interval  for 
error.  ABC  represent  443  222  and  111  observations 
respectively  (after  Carter  &  Robinson  1987) . 


number  of  XBTs  required  in  the  analysis  can  be  determined. 
The  area  of  ocean  under  consideration  in  this  study  is  very 
difficult  to  interpolate  adequately  without  a  large  amount  of 
data  because  of  the  high  degree  of  variability  that  exists 
within  short  spatial  and  temporal  ranges. 

Rather  than  just  looking  at  one  parameter,  such  as  sea 
surface  temperature  or  the  depth  of  the  15  degree  isotherm, 
this  study  will  seek  to  reproduce  the  vertical  temperature 
structure  at  a  given  position  within  the  region  from  the 
surface  to  a  depth  of  8r0  metres.  The  analysis  will  also  give 
an  account  of  the  average  error  involved  in  creating  a 
synthetic  XBT  profile. 

There  are  two  theoretical  strands  that  are  considered;  (1) 
the  theory  of  objective  analysis  and  (2)  the  theory  of 
Empirical  Orthogonal  Functions. 

The  first,  objective  analysis,  describes  a  method  to  take 
a  finite  number  of  data  points,  at  irregular  spatial  or 
temporal  intervals  over  an  area  of  the  ocean's  surface,  and 
interpolate  the  data  in  such  a  manner  that  an  optimal  estimate 
of  a  scalar  value  can  be  obtained  for  any  given  location 
within  the  region. 

The  amount  of  data  to  be  interpolated  per  grid  point  is 
further  reduced  by  exploiting  the  properties  cf  Empirical 
Orthogonal  Functions.  The  use  of  EOFs  allows  a  given  XBT  to  be 
broken  down  into  modes  that  are  constant  for  the  whole  data 
set,  and  into  corresponding  modal  amplitudes  that  are  unique 
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uto  each  particular  XBT.  Then,  for  each  XBT,  the  sum  of  the 
products  of  the  modes  and  the  corresponding  modal  amplitudes 
give  a  complete  representation  of  the  XBT  in  question. 
However,  the  first  few  EOFs  often  explain  the  majority  of  the 
structure  of  the  complete  XBT.  Thus  it  is  possible  to 
approximately  reconstruct  each  XBT  with  a  reduction  in  the 
data.  If,  for  instance,  only  the  first  4  modes  are  considered, 
then  each  XBT  is  represented  by  just  4  unique  numbers,  the 
modal  amplitudes. 

Having  determined  the  unique  modal  amplitudes  for  each 
XBT,  objective  analysis  is  used  to  optimally  interpolate  the 
four  principle  sets  of  numbers  onto  a  regular  grid. 
Multiplying  each  in  turn  by  its  corresponding  mode,  results  in 
a  synthetic  XBT  being  reconstructed  at  each  of  the  grid 
points . 

Having  synthetically  produced  XBTs  at  each  grid  point,  an 
objective  error  analysis  is  used  to  estimate  the  total  error 
variance  of  each  of  the  synthetic  XBTs.  Thereafter,  using  a 
random  generator,  successive  XBTs  are  removed  from  the 
original  data  set.  The  objective  analysis  and  reconstruction 
are  repeated  until  the  error  variance  in  the  synthetic  XBTs 
become  unacceptable. 
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II  THEORY 


A.  INTRODUCTION 

This  chapter  is  divided  into  4  main  sections.  The  first 
section  outlines  the  development  and  theory  of  EOFs  and  gives 
an  account  of  the  use  of  EOFs  in  oceanography.  Similarly, 
section  two  covers  the  background  of  objective  analysis  and  is 
followed  by  a  development  of  the  theory.  The  third  section 
explains  how  the  error  analysis  of  the  modal  amplitudes  is 
used  to  account  for  the  error  in  the  reconstructed  synthetic 
XBT.  The  final  section  outlines  how  all  the  strands  can  be 
brought  together. 

B.  DEVELOPMENT  AND  THECKY  OF  EMPIRICAL  ORTHOGONAL  FUNCTIONS 
1.  Development  of  Empirical  Orthogonal  Functions 

The  theory  of  Principle  Component  Analysis  was  first 
proposed  by  Pearson  (1901)  and  developed  into  a  comprehensive 
theory  by  Hotelling  (1933).  Hotelling's  work  led  Kelly 
(1935)  to  advance  a  model  suitable  for  modern  computer  usage. 
The  theory  was  first  put  into  practice  by  Wrigley  and  Nechus 
(195^)  in  the  field  of  psychology. 

Lorenz  (1956)  outlined  the  theoretical  basis  for  the 
use  of  Principle  Component  analysis  in  meteorology, 
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demonstrating  its  use  as  an  aid  to  efficient  weather 
prediction  and  coining  the  phrase  Empirical  Orthogonal 
Functions  (EOFs)  which  has  become  the  accepted  norm  within  the 
geophysical  sciences.  The  value  of  EOFs  as  a  tool  in 
geophysical  research  is  reflected  in  the  variety  of  uses  to 
which  they  have  been  put.  For  instance,  in  meteorological 
research,  which  requires  working  with  large  data  sets,  the  use 
of  EOFs  have  been  used  to  reduce  the  volume  of  data  that  need 
to  be  interpolated  or  stored. 

Stidd  (1966)  used  EOFs  to  study  climatological  rain 
fall  patterns  within  the  State  of  Nevada.  By  interpolating  EOF 
analysis  between  climate  stations  he  was  able  to  successfully 
reconstruct  the  climate  record  of  a  station  that  had  been 
removed  from  the  initial  analysis.  This  is  similar  to  the 
current  study  in  that  the  temperature  data  at  each  XBT  site, 
like  the  rain  fall  data,  is  represented  by  modal  amplitudes, 
and  the  data  must  be  interpolated  to  additional  locations 
using  objective  analysis. 

EOFs  have  been  featured  highly  in  climatological 
studies  causing  Mitchell  (1966)  to  comment  that  EOFs  may  be  of 
significant  use  as  climatological  indicators.  This  view  is 
strengthened  by  the  work  of  researchers  like  Kutzbach  (1967), 
who  used  EOF  analysis  successfully  to  combine  climatological 
records  of  temperature,  precipitation  and  surface  pressure 
over  the  United  States;  and  Kidson  (1974),  who  used  EOFs  to 
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produce  climatological  indicators  for  both  hemispheres  and  the 
tropics. 

Paegle  and  Haslam  (1982)  used  EOFs  in  the  prediction 
of  the  500  and  850  mb  pressure  heights  over  a  24  hour  period. 
Wallace  and  Dickinson  (1972)  showed  how  EOFs  may  be  applied  to 
time  series  analysis,  reducing  the  data  processing  required 
and  increasing  the  efficiency  for  spectral  modelling  of  the 
atmosphere. 

In  oceanography  the  technique  is  finding  an 
increasingly  wide  variety  of  uses.  For  instance  Kundu  (1975) 
used  EOFs  in  a  time  series  analysis  of  velocity  fields  along 
the  Oregon  coast.  Carnes  et  al  (1990)  have  shown  that  EOFs  can 
be  used  in  conjunction  with  satellite  derived  ocean  dynamic 
heights  to  obtain  a  measure  of  the  ocean's  subsurface  vertical 
temperature  structure.  Oceanographic  models  at  Fleet 
Numerical  Oceanographic  Center  (FNOC) ,  such  as  the  Optimal 
Thermal  Interpolation  System  (OTIS) ,  employ  EOFs  to 
effectively  represent  ocean  thermal  climotologies 
(Tunnicliffe  and  Cummings  1991) .  Similarly,  the  Navy/NOAA 
Oceanographic  Data  Distribution  System  (NODDS)  includes  the 
use  of  EOF  techniques  to  compress  large  volumes  of  data, 
enabling  distant  users  either  ashore  or  at  sea  to  receive  by 
telephone  link  sophisticated  real  time  ocean  and 
meteorological  information  using  a  desk  top  PC. 
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2.  Theory  of  Empirical  Orthogonal  Functions's 


The  above  examples  show  the  versatility  and  value  of 
EOFs  as  an  effective  tool  within  the  fields  of  meteorology  and 
oceanography.  Set  out  below  is  a  development  of  the  basic 
theory.  The  approach  outlined  considers  the  work  of  Lorenz 
(1956),  who  first  described  the  use  of  EOFs  in  geophysical 
research,  Harman  (1976),  who  formally  derives  the  general 
theory,  and  Dunteman  (1989),  whose  clarity  and  examples  gave 
considerable  insight  into  the  technique. 

The  object  of  Principle  Component  Analysis  is  to  take 
a  large  body  of  data  and  empirically  reduce  it.  The  model 
assumes  a  linear  set  of  numbers  such  that  a  linear  combination 
of  these  components  leads  to  a  complete  representation  of  the 
original  data  set. 

Mathematically  the  method  assumes  that, 


^i=<7i>ri+<3r2>r2+ . Qjyj 

equation  1 

where  (  j  =  1,2, . n)  and  each  of  the  observed  variable  Pl 

is  described  linearly  in  terms  of  n  orthogonal  components 
y1,yz,.,.yn.  The  power  of  this  approach  being  that  only  a  few 

of  the  components  need  to  be  retained  in  order  to  retain  the 
majority  of  the  total  variance. 
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The  coefficients  g^ ,  are  referred  to  as  the 

"loadings,  "scores"  or  "weightings"  and  in  geophysics  as 
"modal  amplitudes".  Each  modal  amplitude  is  multiplied  by  its 
corresponding  principle  component,  with  the  sum  being  equal  to 
the  value  of  the  original  variable.  The  problem  is  to  find 
suitable  values  of  q  and  y  to  be  able  to  represent  the 
variable  pi  in  question.  This  is  efficiently  achieved  by 

expressing  equation  1  in  matrix  form  as. 


P=QY 


equation  2 


where  P  is  an  m  by  n  matrix  of  scalar  variables  whose  columns 
represent  the  vector  Pi  ,  Q  is  an  m  by  n  matrix  of  modal 

amplitudes  and  Y  represents  n  column  vectors  each  with  m  rows. 

Consider  the  situation  where  m  elements  (Pii=l...m) 

have  been  measured  at  n  different  locations.  For  this  study, 
80  isotherm  depth  measurements  (m)  made  at  each  of  156 
locations  (n) . 

Let 


a=pV 

equation  3 
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where 


P*=Pi-P 

equation  4 

the  difference  of  a  value  Pi  from  the  mean  value  T.  P *“  is 

the  transpose  of  P*  and  A  represents  the  covariance  matrix 
formed  by  the  dot  product  of  P*‘  and  P* .  The  covariance  matrix 
A  is  normalized  to  form  the  correlation  matrix  A, 

A  =  <P*/XP*> 

equation  5 

and  the  symbols  <>  denote  an  ensemble  averaging  of  the 
variance  from  each  data  point.  The  matrix  A  is  also  known  as 
the  projection  matrix. 

From  the  theory  of  matrix  algebra  (Harman  1976)  a 
general  matrix  O  can  be  expressed  in  terms  of  its  eigenvectors* 
and  eigenvalues  X  such  that, 


Oe-Xe 


equation  6 


12 


where  09  is  regarded  as  a  transformation  of  0  with  X  as  the 
constant  of  proportionality.  Each  root  Xi  has  a  non  zero 

solution  et  and  the  m  roots  Xx,X2, . Xm  lead  to  n  values 

ex,e2.....en  such  that  equation  6  may  be  written  as, 


O  (  ,  &2  ....  '  ^2  ^2 


en 


or  in  matrix  form 


gb=bA 


equation  7 


where  A =diag(X1,  Xz.  .  .  .  Xn)  . 


equation  8 


Inserting  matrix  A  from  equation  5  into  the  general 
equation  8  gives, 


AY=XY 


equation  9 


The  vectors  yi  are  linearly  independent  such  that  the 
determinants  are  non  zero  and  Y  has  a  unique  inverse  3TX, 


y^ay=A. 


equation  10 
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Since  A  is  a  correlation  matrix  and  is  symmetric 


A=Af 

equation  11 

(Y~x)'  are  characteristic  values  with  Y  being  orthogonal  with 
the  property  that, 


YY*=I 


equation  12 


the  identity  matrix  or. 


Yjl=Y' 


equation  13 


giving  that  equation  14  can  be  written  as, 

Y*AY=  A 


equation  14 


Equation  14  states  that  the  symmetric  matrix  A  may 
be  diagonalized  by  means  of  the  orthogonal  transformations Y 
and  that  the  elements  of  A  and  Y  are  real  with  Y  being  made 
up  of  n  characteristic  linearly  independent  equations. 
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From  equation  14  diagonally  decomposing  A  gives 
values  for  matrix  Y,  and  A.  Now  knowing  the  values  of  P  and 
Y ,  Q,  the  matrix  of  coefficients  or  modal  amplitudes,  can  be 
determined  from  equation  2.  Having  obtained  values  for  Q  andlT 
equation  2  states  that  the  value  of  P  can  be  exactly 
determined  and  that  it  equals  the  matrix  product  of  Q  and  Y. 

From  equations  2, 13, and  14  it  follows  (Paegle  and 
Haslam  1982)  that  the  total  variance  is  given  by  the  sum  of 
the  eigenvalues. 


equation  15 

and  that  each  eigenvalue  Xi  gives  the  contribution  of  each 
eigenvector  Ys  to  the  total  variance  of  P. 

When  the  eigenvalues  are  arranged  in  descending  order 
the  variance  represented  by  each  mode  or  eigenvector  decreases 
dramatically  as  the  number  of  the  eigenvalues  are  increased. 
A  realistic  estimate  of  the  original  data  P  can  thus  be 
achieved  by  using  only  the  first  few  modes  Y  and  the 
corresponding  modal  amplitudes  Q. 
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P=QY/. 


equation  16 

The  use  of  a  limited  number  of  modes  reduces  the 
quantity  of  data  that  has  to  be  stored  and  processed.  In 
addition  the  variability  in  the  higher  modes  is  likely  to 
represent  noise  in  the  original  signal.  Thus,  by  removing  the 
higher  modes  a  "cleaner  profile"  is  obtained.  Preisendorfer 
et  al  (1981)  suggest  that  modes  which  can  not  be  distinguished 
from  randomly  generated  data  should  be  removed.  Dunteman 
(1989)  suggests  that  all  modes  for  which  the  eigenvalue  is 
less  than  one  should  be  removed.  Dunteman' s  approach  is  used 
within  this  study. 

C.  DEVELOPMENT  AND  THEORY  OP  OBJECTIVE  ANALYSIS 
1.  Development  of  Objective  Analysis 

Objective  analysis  is  a  technique  that  will  produce 
an  optimal  estimate  of  some  quantity  at  a  given  location  by 
the  interpolation  of  irregularly  spaced  data  points.  The 
method  is  based  upon  the  Gauss  -  Markov  theory . 

Objective  analysis  was  first  used  in  meteorology  by 
Gandin  (1965)  who  used  the  technique  to  analyze  atmosphe:  o 
pressure  and  windfields.  The  technique  was  introduced  for 
oceanographic  use  by  Bretherton  et  al  (1976) ,  who  demonstrated 
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its  value  in  determining  optimal  temperature,  velocity  and 
streamline  maps.  The  technique  was  applied  by  Freeland  and 
Gould  (1976)  to  data  taken  during  POLYMODE  and  successfully 
produced  stream  function  maps  of  the  North  West  Atlantic. 

Carter  (1983)  extended  the  use  of  objective  analysis 
by  considering  distance  variations  separately  in  the  X  and  Y 
directions  and  a  temporal  component,  thereby  allowing 
observations  made  at  different  places  and  at  different  times 
to  be  mapped.  In  addition,  the  theory  allows  an  explicit 
statement  to  be  made  about  the  error  in  the  determination  of 
an  interpolated  value  at  a  given  location.  Because  of  the 
introduction  of  a  temporal  component,  Carter's  method  also 
enables  maps  of  the  quantity  to  be  predicted  for  a  future 
time. 

Objective  analysis  is  now  widely  used  in  oceanography. 
For  instance.  Watts  et  al(l989)  used  objective  analysis  to 
model  the  depth  of  the  12  degree  Celsius  isotherm  from 
inverted  echo  sounder  observations  taken  in  the  vicinity  of 
the  Gulf  Stream.  Objective  analysis  is  a  standard 
interpolation  tool  that  is  extensively  used  for  computer 
aided  numerical  prediction  in  both  meteorology  and 
oceanography  (see  Clancy  1989) . 

2.  Theory  of  Objective  Analysis 

The  derivation  outlined  below,  after  Carter  (1983), 
forms  a  statement  of  the  Gauss  Markov  theory  for  determining 
a  least  squares  optimal  value. 
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The  statistical  model  for  objective  analysis  assumes 
a  stationary  homogeneous  field.  Let  @r  be  a  measurement  of 

some  quantity  and  let  the  error  in  the  measurement  be  er . 
Then, 


0r=9r+er 


equation  17 


where  0r  is  the  true  value.  It  is  assumed  that  observation 
error  is  uncorrelated  with  the  true  field  such  that, 

R(ezBs)  =0 

equation  18 


where  J?(er8x)  represents  the  correlation  between  the  error  e 

at  position  r  and  the  measured  field  at  some  other  locations. 

It  is  also  assumed  that  the  correlation  between 
observation  errors  at  two  locations  is  zero, 


R<.eres)  =e261 


equation  19 


where  R{eres)  represents  the  correlation  between  er  and  e,,e2 
is  the  error  variance,  and  6  is  the  Krondiker  delta  having 
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a  value  of  one  when  r  equals  s  and  the  value  of  zero 
otherwise. 

Objective  analysis  seeks  to  find  the  optimal  value  of 
a  given  quantity  X  at  an  arbitrary  location.  The  optimal 
estimate  of  the  value  at  the  grid  location  is  designated 
In  matrix  form  the  estimate  at  the  grid  points  is  given  as  a 
linear  combination  of  the  values  of  the  data  measured  at  a 
variety  of  locations  r  such  that, 

X-A0r 

equation  20 

where  0^  is  the  value  of  the  quantity  measured  at  position  r 
throughout  the  region.  For  example  0r  could  represent  sea 

surface  temperature  measurements  taken  at  various  irregularly 
spaced  positions  within  a  given  region.  Whereas  X  represents 
true  values,  the  value  X  is  the  estimate  that  is  determined 
at  the  grid  points  by  interpolating  the  values  of  0  onto  the 
grid  by  the  use  of  linear  combinations  of  0^.  using  the  matrix 
A. 

In  order  to  determine  the  estimates  at  the  chosen 
gr..d  points  it  is  first  necessary  to  ascertain  values  for  the 
elements  in  matrix  A .  This  is  done  in  such  a  manner  as  to  give 
the  optimal  estimate  of  2.  Throughout  the  derivation  X  andJ? 
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are  referred  to  as  if  they  were  known,  whereas  in  fact  they 
are  the  quantities  ultimately  that  are  to  be  determined. 

Initially  it  is  the  value  of  A  that  is  sought  such 
that  it  minimizes  the  error  by  a  least  squares  fit  between  the 
true  value  of  X  at  the  grid  points  and  the  estimate  X. 
Firstly  let 


CM=E[X 6'] 


equation  21 


where  C*  is  the  correlation  matrix  found  by  comparing  the 

value  of  the  quantity  at  the  required  grid  point  locations 
compared  with  those  at  the  given  data  sites. 

Cx~E[XX/] 


equation  22 


where  Cx  is  the  correlation  between  the  value  at  any  required 

grid  point  location  compared  to  the  value  at  any  other 
required  grid  point  location. 


q,=i?[e0'] 


equation  23 


where  C*  is  the  correlation  between  the  values  at  any  two  data 
point  sites. 
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Then  to  obtain  the  optimal  interpolation  the  value  of 
the  error  C9  is  minimized  such  that 

Ct=E[ee']  =E[(J?-X)  (Jt-X)'] 

equation  24 

where  C0  represents  the  correlation  between  the  mean  square 

variance  of  the  estimated  values  compared  to  the  actual 
values.  Substituting  equation  20  into  equation  23  and 
expanding  gives, 

C0=E[  (A9-0)  <A0-O)') 

equation  25 


and, 

C.=AC0A/-CrfA/-AC^+ Cx . 

equation  26 

This  expression  can  be  simplified  by  using  a  matrix  identity 
and  noting  that  , 


C#=  (A^CjjQj  )  C0  ( A~ )l~Ca^C^  Cjj + Cx , 

equation  27 
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Since  the  matrices  Cg  and  C£  are  nonnegative  definite,  then 


the  error  matrix  is  minimized  when, 

A-C^C^O 

equation  28 

giving, 

equation  29 

The  value  of  the  error  matrix  C#  can  be  written  explicitly  as, 

— Ox~ C^q Cq 

equation  30 

From  equation  29  and  substituting  for  A  in  equation  20,  the 
estimate  of  the  value  of  the  quantity  at  the  grid  points  is 
given  by, 
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equation  31 

and  the  error  in  these  estimates  is  given  by  equation  30. 

Thus,  providing  the  correlation  matrices  andCg 

can  be  determined  a  value  for  2,  the  estimate  of  the  value  at 

any  given  grid  location  can  be  obtained  from  a  knowledge  of  8, 
the  value  at  any  given  location.  Equations  29  and  30  are  a 
statement  of  the  Gauss  Markov  theory. 

3.  The  correlation  function 

The  correlation  matrix  C^,  a  measure  of  the 

correlation  between  the  values  at  each  of  the  data  sites 
compared  to  the  values  at  each  of  the  grid  points,  is  unknown. 
Similarly,  the  correlation  matrix  Cx,  the  correlation  between 

successive  grid  point  values,  is  also  unknown.  The  only 
correlation  that  is  available  is  Cg,  the  correlation  between 

data  values  at  the  irregularly  sampled  locations. 

However,  the  determination  of  Gg  is  not  straight 

forward.  In  order  to  determine  the  correlation  between  two 
points  it  is  necessary  to  have  made  several  readings  at  each 
location,  whereas  in  this  study  only  one  reading  at  a  given 
location  is  available.  This  problem  is  overcome  by  assuming 
that  the  correlation  between  any  two  points  is  a  function  of 
distance. 
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The  correlation  matrix  is  formed  by  computing  the 

distance  between  each  data  point  and  every  other  data  point. 
The  data  pairs  are  grouped  into  distance  bins,  and  the 
correlation  between  distance  bins  is  then  determined  using  the 
expression, 


6  E  (9,-8g  <6, -Bp 

"  <E  <9r-8PJE  <a#-sp2!I/J 

equation  32 

where  0r  and  8S  are  data  values  at  two  points  r  and  s,  and  Wk 

is  the  mean  of  the  values  for  distance  bin  k.  Once  the 
correlation  function  has  been  determined  for  the  data  points 
within  the  region  the  results  are  applied  to  the  two  unknown 
matrices  and  Cx.  Simply  knowing  the  distance  between  a 

grid  point  and  a  data  point  or  between  two  particular  grid 
points  is  sufficient  information  to  enable  the  corresponding 
correlation  between  the  two  points  to  be  computed. 
Unfortunately  there  is  one  more  slight  complication,  in  that 
the  two  matrices  have  to  be,  by  definition,  positive  definite 
for  equation  31  to  be  valid.  This  means  that  an  estimate  of 
the  correlation  between  two  successive  points  can  not  be 
achieved  from  a  database  simply  by  interpolating  between  two 
adjacent  distance  bins,  because  the  approximation  may  not  be 
positive  definite,  in  order  to  ensure  that  the  two  matrices 
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are  positive  definite  it  is  necessary  to  fit  a  function  to  the 
distance  correlation  database. 

The  function  that  is  normally  fitted  to  the  curve 
(Carter  and  Robinson  1987)  takes  the  form, 

CIS=(  i-{r/a)2)e-,r/w2 

equation  33 

where  a  and  b  are  the  unknowns  to  be  determined,  r  the 
distance  between  any  two  data  points  r  and  s,  and  Cia  the 

correlation  between  them.  The  values  of  a  and  b  are 
determined  iteratively  by  minimizing  the  error  between  the 
original  correlations  CD  as  given  in  the  database  outlined 

above  and  Cr  a .  where  the  error  is  given  by, 


*=(CIS-Cn)\ 

equation  34 

The  correlation  matrices  of  C&  and  Cx  can  now  be 

determined  from  the  function  described  in  equation  33  and  the 
objective  analysis  can  be  undertaken. 
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D.  INTERPOLATED  ERROR 


Having  obtained  an  estimate  for  the  value  of  the  first 
four  modal  amplitudes  at  each  of  the  grid  positions,  it  then 
remains  to  use  the  theory  of  EOFs  to  reconstruct  a  synthetic 
XBT  at  each  of  these  positions.  This  is  simply  achieved  by 
multiplying  each  modal  amplitude  by  its  corresponding 
eigenvector  and  adding  the  four  resulting  vectors  together,  as 
per  equation  16.  However  some  of  the  estimated  modal 
amplitudes  used  contain  error.  The  error  variance  of  each 
modal  amplitude  is  specified  by  equation  30,  and  must  be  taken 
into  account  in  reconstructing  a  synthetic  XBT  at  a  grid 
position. 

Consider  a  modal  amplitude  at  a  particular  grid  location 
Q±  having  an  error  variance  6q,  and  assuming  the  synthetic 

XBT  at  that  position  is  going  to  be  reconstructed  using  i 
EOFs,  then  the  error  variance  in  the  synthetic  XBT  el  can  be 

shown  (Carter  1983)  to  be  given  by, 


equation  35 

The  error  variance  from  this  reconstruction  is  then  mapped 
to  give  a  pictorial  image  of  areas  within  the  region  that 
have  high  and  low  error  variances. 
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The  error  variance  is  a  measure  of  the  confidence  of  a 
given  reconstruction.  Figure  2  shows  an  example  of  an  error 
variance  map.  Low  confidence  is  indicated  when  the  values 
approach  one.  This  map  is  the  combination  of  the  individual 
modal  amplitude  error  maps  shown  in  Figure  3  using  equation 
35.  The  figure  also  shows  where  each  XBT  cast  was  taken.  As 
would  be  expected,  the  lowest  error  variance  (highest 
confidence)  occur  in  areas  that  have  a  high  number  of 
samples,  with  the  error  variance  (lowest  confidence)  being 
largest  where  there  are  no  or  few  samples. 


E.  APPLICATION  OP  THEORY  TO  CURRENT  STUDY 

All  the  elements  of  the  theory  can  now  be  put  together  to 
analyze  the  area  under  investigation.  Firstly  the  original 
XBTs  will  be  converted  into  a  correlation  matrix,  where  one 
depth  is  compared  to  another  and  the  whole  data  set  ensemble 
averaged  to  give  the  projection  matrix.  This  matrix  will  then 
be  decomposed  to  find  the  significant  eigenvectors,  noting  the 
value  of  the  corresponding  eigenvalues.  The  most  significant 
eigenvectors  or  modes  will  be  selected,  and  for  each  XBT 
within  the  set  the  corresponding  modal  amplitudes  will  be 
determined. 

Once  the  modal  amplitudes  have  been  found,  a  correlation 
matrix  as  a  function  of  distance  can  be  constructed.  From  this 
an  appropriate  function  will  be  fitted  and  the  correlation 
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matrices  Cx  and  CjO  determined.  The  modal  amplitudes  can 

then  be  optimally  interpolated  onto  a  grid.  The  process  is 
repeated  for  the  second,  third,  and  fourth  modal  amplitudes. 

Synthetic  XBTs  can  then  be  reconstructed  at  each  grid 
point  using  the  interpolated  modal  amplitudes  and  a  measure  of 
the  error  variance  in  each  XBT  can  be  determined  from  the 
error  matrices  generated  by  the  objective  analysis. 


Figure  2  Reconstructed  error  variance  map  using  all  156 
XBT's.  The  map  is  produced  using  equation  37  (effectively 
combining  the  four  maps  from  Figure  2.  The  contours  are  at 
0.1  spacing.  The  central  contour  represents  0.1  (or  10%) 
error  variance. 
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Figure  3  Error  variance  maps  for  the  first  four  modal 
amplitudes.  Top  left  shows  error  variance  for  first,  top 
right  for  second,  bottom  left  and  right  show  third  and 
fourth  respectively.  At  0.1  intervals.  Inner  most  contour 
representing  0.1  or  10%  error  variance. 


Ill  DATA 


A.  THE  MEANDER  EXPERIMENT 

The  data  for  this  study  consists  of  156  XBTs  taken  within 
the  region  of  a  Gulf  stream  meander  sampled  during  the  period 
September  17th  to  October  13th  1988.  The  original  data  was 
collected  as  part  of  a  much  wider  experiment  that  involved  two 
cruises,  one  in  the  autumn  of  1988  and  the  second  in  the 
spring  of  1989.  The  first  cruise  sampled  an  anticyclonic 
meander  crest  (EN  185)  where  the  path  of  the  current  is  convex 
to  the  North.  The  second  cruise  collected  data  from  a  cyclonic 
meander  trough  (EN  194)  where  the  flow  of  the  current  is 
convex  to  the  South.  The  objective  of  the  two  cruises  was  to 
investigate  the  time  dependent  kinematics  and  dynamical 
structures  of  Gulf  Stream  meanders.  The  Gulf  Stream  meander 
was  sampled  with  a  variety  of  instruments,  and  density  and 
velocity  fields  were  computed  to  enable  fluxes  of  mass, 
momentum  and  vorticity  to  be  determined  as  the  meander 
progressed  in  space  and  time. 


B.  THE  XBTs 

The  following  technical  details  of  the  XBTs  are  taken  from 
the  initial  cruise  report  (Hummon  1991) . 
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The  XBTs  used  in  the  survey  were  Sippican  T7  probes  which 
have  a  nominal  depth  rating  of  7  60  metres.  The  XBTs  were 
launched  from  a  fixed  stern  deck  launcher  with  a  BathySystemn 
810  XBT  deck  unit.  The  data  was  stored  on  a  HP-85B  computer 
equipped  with  an  HP9121D  disk  drive.  The  software  was  supplied 
by  BathySystems  but  was  substantially  modified  to  allow 
simpler  and  faster  processing.  The  raw  data  was  recorded  in 
volts  versus  descent  time,  The  data  were  transferred  to  a 
MassComp  computer  and  each  profile  was  converted  into 
temperature  versus  depth  measurements  and  stored  onto  disk  or 
magnetic  tape. 

The  resolution  of  the  data  is  0.65  metres  with  a  0.1  metre 
precision.  The  stated  accuracy  of  the  depth  measurement  is 
five  metres  or  2%  of  the  depth,  whichever  is  greater. 
Temperature  data  is  stored  to  within  0.001  degree  Celsius 
with  measurement  accuracy  to  within  0.15  degrees  Celsius. 

The  data  was  edited  to  remove  the  first  three  measurements 
corresponding  to  depths  less  than  two  metres.  Readings  taken 
at  depths  greater  than  810  metres,  outside  of  the  stated 
operating  range  of  the  probes,  were  also  removed.  Spikes,  bad 
data  and  wire  breaks  in  individual  profiles  were  deleted  by 
hand  on  the  MassComp  computer.  The  full  set  of  XBT  casts  is 
shown  in  Figures  4-17.  The  geographic  distribution  of  the 
casts  is  shown  in  Figure  2,  with  location  values  being  given 
in  the  log  shown  in  Appendix  A. 
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Figure  4  -  17  show  all  the  XBT  casts  taken  during  the 
Anatomy  of  a  Meander  experiment. 
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Figure  6  XBTs  25-36 


Figure  7  XBTs  37  -  48 


Figure  8  XBTs  49  -60 
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Figure  10  XBTs  73  -84 
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Figure  n  XBTs  85-96 
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Figure  12  XBTs  97  -  108 
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Figure  14  XBTs  121  -  132 
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Figure  17  XBTs  157  -  168 
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IV  METHODS 


Part  A  of  this  chapter  describes  how  the  original 
projection  matrix  was  computed  and  how  the  eigenvectors  and 
eigenvalues  were  determined.  Part  B  describes  how  the 
Objective  Analysis  was  implemented  and  how  the  synthetic  XBTs 
were  reconstructed.  Finally  part  C  describes  how  the  data  set 
was  reduced  to  find  the  minimum  number  of  XBT  sites  that  were 
required  in  the  case  of  this  particular  Gulf  Stream  meander. 


A.  DEPTH  CORRELATION  MATRIX 
1.  The  matrix 

To  overcome  initial  data  analysis  problems  all  XBTs 
less  than  800  metres  were  removed  from  the  data  set.  This  left 
a  total  of  156  useable  XBTs  for  further  analysis. 

The  vertical  correlation  matrix  was  formed  using 

FORTRAN  program  LOADBATHYS  (Appendix  IB)  and  subroutine  REDATA 

* 

(Appendix  2B) .  The  subroutine  interpolates  temperature  values 
from  each  XBT  at  10  metres  intervals  commencing  with  a  depth 
of  five  metres.  The  vertical  correlation  matrix  was  computed 
in  the  main  program  by  comparing  the  temperature  at  one  depth 
with  that  at  another  depth.  This  process  ensemble  averaged 
over  all  156  XBTs  using  equation  36. 
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*=£  a>  1/3 


equation  36 

where  A  is  the  80  by  80  projection  matrix  formed  by  comparing 
the  temperature  at  all  80  depths  with  each  other  and  8  is  the 
temperature  at  depths  i  and  j ,  with  the  overbar  representing 
the  mean  temperature  for  that  depth  i  or  j .  The  projection 
matrix  is  visualized  in  Figure  18. 


temperature  between  depths. 


2.  Eigenvectors  and  values 


The  correlation  matrix,  A,  was  decomposed  to  find  its 
eigenvalues  and  eigenvectors  using  equation  9. 

Figure  19  shows  the  first  six  eigenvalues  and  the 
associated  variance  for  the  first  four  modes.  Each  eigenvalue 
is  proportional  to  the  variance  contributed  by  its 
corresponding  eigenvector  (Harman  1976) .  The  first  eigenvector 
accounts  for  over  75%  of  the  variance  of  the  correlation 
matrix,  the  second  eigenvector  is  responsible  for  15%,  the 
third  for  5.1%, and  the  fourth  for  1.7%.  The  cumulative  percent 
variance  explained  by  the  first  four  eigenvectors  is  over  98% 
of  the  total  variance  of  the  projection  (correlation)  matrix. 
Thus,  instead  of  using  80  eigenvectors  to  describe  the 
variance  in  the  correlation  matrix  A,  it  is  possible,  using 
the  criteria  discussed  by  Dunteman  (1989),  to  describe  the 
matrix  sufficiently  with  only  four,  with  a  minimal  loss  in 
information,  thereby  saving  considerably  on  data  storage  and 
processing  requirements  and  suppressing  the  noise  contained 
within  the  higher  modes.  The  modal  amplitudes  for  each  XBT 

i 

were  calculated  using  equation  16  in  a  MATLAB  subroutine. 
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Figure  19  The  first  6  eigenvalues.  Percentage  of  variance 
is  shown  for  first  4  modes. 


B.  OBJECTIVE  ANALYSIS 

Each  XBT,  after  the  application  of  the  EOF  decomposition, 
is  represented  by  four  modal  amplitudes.  The  problem  is  to 
interpolate  ;he  modal  amplitudes  to  arbitrary  positions 
within  the  analysis  region  using  objective  analysis.  It  was 
decided  to  compute  the  interpolated  modal  amplitudes  at 
regular  intervals  using  a  half  degree  spacing  in  both 
longitude  and  latitude,  over  a  grid  extending  from  37  to  40 
degrees  North  and  64  to  72  degrees  West.  The  length  scale 
between  grid  points  of  approximately  50  km  was  chosen  because 


9 


it  is  comparable  to  the  Rossby  radius  of  deformation  at  this 
latitude. 

The  estimate,  2,  of  each  amplitude  at  each  grid  location 
using  equation  31  was  computed.  The  correlations  are  assumed 
to  depend  solely  upon  the  distance  between  observations  and 
similarly  between  observations  and  grid  points. 

1.  Determination  of  special  correlation  matrices 

The  distance  between  XBT  sites  was  calculated  and 
grouped  into  bins.  Several  bin  intervals  were  considered  with 
the  object  being  to  find  an  interval  that  gave  a  reasonable 
number  of  data  pairs  per  bin,  allowing  an  unbiased  measure  of 
correlation  by  distance  to  be  determined.  This  was  achieved 
using  the  program  DEEPCOR  and  the  subroutine  CALC  described  in 
Appendix  3B. 

The  number  of  data  pairs  for  the  three  intervals  used 
are  shown  at  Table  I.  The  25  km  interval  gave  sufficient 
data  pairs  for  each  bin  out  to  a  distance  of  200  km,  and 
allowed  eight  spatial  correlation  estimates  to  be  made.  The 
results  are  shown  in  Figures  20-23. 

4 

The  correlation  function  described  in  equation  33  was 
used  to  model  the  correlation  estimates  shown  in  Figures  20- 
23.  The  parameter  a  is  equal  to  the  distance  at  which  the 
correlation  falls  to  zero,  and  is  given  as  the  point  where  the 
curve  in  Figures  20-23  crosses  the  X  axis.  The  value  of  b  is 
the  value  of  the  distance  when  the  correlation  equals  the  e 
folding  distance  (e’1).  The  value  of  the  coefficients  a  and  b 


50 


were  found  iteratively  using  the  program  FUNCTION  given  in 


Appendix  4B.  In  this  program,  the  square  error  (,), 


«=(Cr-Cn) 


equation  37 


between  the  original  data  points,  r  ,  in  Figures  20-23,  and 

the  iterated  values  q  ,  calculated  using  equation  33,  is 

minimized.  The  iteration  sequence  is  intialized  with  values 
of  a  and  b  from  visually  inspecting  Figures  20-23. 

In  order  to  determine  whether  each  incremented  value 
of  a  and  b  should  be  larger  or  smaller  than  the  initial 
value,  equation  33  was  differentiated  with  respect  to  a  and 
with  respect  to  b.  The  analytical  solution  was  used  to 
increment  a  and  b  in  such  a  way  that  the  mean  square  error  was 
reduced  with  each  iteration.  The  iteration  was  repeated  until 
the  error  had  reduced  to  0.05.  The  final  values  of  the 
parameters  a  and  b  are  shown  in  Table  II. 

It  is  assumed  that  the  distance  correlation  function 
determined  above  for  ^  will  also  be  applicable  in  the 

observation  to  grid  point  correlation  matrix  q  . 

The  objective  analysis  FORTRAN  source  programs  are 
provided  in  Appendix  5B,  6B  and  7B  for  reference.  The  first 
guess  for  each  analysis  is  taken  as  the  local  weighted  average 
of  the  modal  amplitudes.  Output  from  the  objective  analysis 
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consists  of  contour  maps  of  the  first  four  modal  amplitudes, 
and  analysis  error  of  the  interpolated  amplitudes. 

Contour  plots  of  the  first  four  modal  amplitudes  are 
shown  in  Figures  24-27  and  their  associated  error  maps  in 
Figures  28-31. 

2.  The  reconstruction 

Using  the  modal  amplitudes  calculated  by  the 
objective  analysis,  synthetic  XBTs  were  reconstructed  at 
each  of  the  grid  points  using  Equation  16.  However,  the 
error  in  the  XBT  reconstruction  is  dependent  on  the  position 
of  the  reconstruction.  Synthetic  XBTs  produced  in  areas  with 
high  concentrations  of  observation  stations  are  expected  to 
suffer  less  error  in  reconstruction  than  synthetic  XBTs 
produced  in  areas  with  sparsely  populated  data.  The  error 
variance  in  each  XBT  was  calculated  using  equation  34  and  the 
resulting  error  variance  map  is  shown  in  Figure  32. 

C.  REDUCING  THE  NUMBER  OF  XBTs 

Of  ultimate  interest  is  the  size  of  the  error  variance  in 
XBTs  reconstructed  within  the  analysis  area.  From  the 
associated  error  variance  map  it  is  possible  to  assess,  for 
any  given  position,  the  value  of  reconstructing  and  using  a 
synthetic  XBT  at  that  point. 

It  was  decided  that  for  a  reconstructed  synthetic  XBT, 
less  than  30%  error  could  be  of  use.  The  area  inside  the  30% 
contour  of  Figure  32  was  noted.  Successive  XBTs  were  removed 
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and  the  objective  analysis  repeated  until  the  30%  contour 
became  the  central  or  first  contour.  This  meant  that  the  area 
that  was  now  enclosed  represented  error  variances  greater 
than  20%  but  less  than  30%.  The  number  of  XBTs  remaining  was 
noted. 

The  original  XBTs  were  numbered  sequentially  and  the 
FORTRAN  program  RANDUM  was  used  to  place  these  numbers  in 
random  order.  On  commencing  the  objective  analysis  suite  of 
programs,  subroutine  REDUCE  permitted  the  number  of  XBTs  to 
be  used  in  the  objective  analysis  to  be  varied. 

Table  I  NUMBER  OF  DATA  POINT  PAIRS  PER  BIN  FOR  THREE 
DIFFERENT  BIN  SIZES. 


12  Km  bin  size 

12 

24 

36 

48 

60 

72 

84 

96 

108 

120 

81 

104 

114 

142 

184 

190 

214 

226 

210 

168 

j  25  Km  bin  size 

25 

50 

75 

100 

125 

150 

175 

200 

225 

250 

195 

266 

412 

452 

370 

310 

276 

224 

90 

68 

|  50  Km  bin  siz< 

a 

50 

100 

150 

200 

250 

300 

350 

400 

450 

500 

461 

864 

680 

500 

158 

58 

52 

62 

92 

106 
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Figure  20  Correlation  between  data  point  pairs  using  a  25 
km  distance  bin  for  the  first  modal  amplitude. 


Figure  21  Correlation  between  data  point  pairs  using  a  25 
km  bin  for  the  second  modal  amplitude. 


100  ISO  200 

distance  25  km  bins 


Figure 
km  bin 


22  Correlation  between  data  point  pairs  using  a  25 
for  the  third  modal  amplitude. 


100  ISO  200 

distance  25  tan  bins 


Figure  23  Correlation  between  data  point  pairs  using  a  25 
km  bin  for  the  fourth  modal  amplitude. 
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Table  II  PARAMETER  VALUE  a  AND  b  FOR  EACH  OF  THE  MODAL 
AMPITUDES . 


Figure  24  Contour  map  of  the  first  modal  amplitude 


Figure  25  Contour  map  of  the  second  modal  amplitude. 
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Figure  27  Contour  map  of  the  fourth  modal  amplitude. 


Figure  28  Error  variance  map  as  a  result  of  interpolating 
the  first  modal  amplitudes.  0.1  contour  intervals. 


Figure  29  Error  variance  map  as  a  result  of  interpolating 
second  modal  amplitudes.  0.1  contour  intervals. 


Figure  30  Error  variance  map  as  a  result  interpolating 
third  modal  amplitudes,  o.l  contour  intervals 


Figure  32  Error  variance  map,  produced  by  using  equation 
37,  effectively  the  combination  of  Figures  28-31.  The  * 
indicates  the  original  XBT  sites.  The  contours  represent 
the  amount  of  confidence  that  can  be  placed  in  a 
reconstruction  at  any  location  within  the  area. 
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V  INITIAL  ANALYSIS 


A.  RECONSTRUCTION  ALONG  A  LINE  OF  LATITUDE 

To  get  a  feel  for  how  good  or  bad  the  reconstructions 
appeared,  a  line  of  synthetic  XBTs  along  37.5  N  from  72  West 
to  69.0  West  were  reconstructed  at  1/2  degree  intervals  and 
are  shown  in  Figures  33-36.  A  group  of  real  XBTs,  taken  along 
the  proximity  of  this  line  are  shown  in  Figures  37-  44.  The 
positions  of  the  real  XBTs  are  readily  apparent  by  consulting 
Figure  45  which  indicates  the  position  of  the  XBTs  used  in 
this  analysis. 

Although  the  two  groups  of  diagrams  show  a  general 
similarity  in  shape  there  are  enough  differences  to  cause 
concern.  Firstly,  the  synthetic  XBTs  all  tend  to  exhibit  an 
temperature  minimum  at  about  200  metres  that  is  more 
exaggerated  than  in  the  real  XBTs  surrounding  this  line  of 
latitude,  secondly,  the  synthetic  XBTs  also  show  a  strong 
negative  temperature  gradient  within  the  first  30  to  80  metres 
that  again  is  not  apparent  in  most  of  the  real  XBTs,  which 
for  the  most  part  are  isothermal  or  exhibit  only  a  slightly 
negative  temperature  gradient  over  the  same  depth  range. 
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B.  RECONSTRUCTION  OF  ONE  XBT 


To  pursue  these  discrepancies  further  XBT  7105A  was 
selected  for  closer  study.  This  particular  XBT  was  chosen 
because  its  position  is  at  the  same  location  as  a 
synthetically  produced  XBT.  Thus,  it  would  be  expected  that 
the  profiles  of  the  real  and  the  synthetic  XBTs  should  show  a 
very  high  degree  of  similarity.  The  real  XBT  7105A,  is  shown 
in  Figure  42  and  the  synthetic  XBT  in  Figure  34.  Again,  the 
synthetic  profile  exhibits  a  temperature  minimum  at  two 
hundred  metres  and  a  negative  gradient  in  the  surface  layer, 
both  features  being  less  pronounced  in  the  observed  profiles. 

As  a  check  to  ensure  that  the  EOF  decomposition  had  been 
performed  correctly  it  was  decided  to  reconstruct  7105A  using 
all  80  modes.  The  results  of  this  are  shown  in  Figure  46  where 
comparison  with  the  original  and  the  synthetic  XBT  using  four 
modes  can  be  made  directly.  Figure  47  shows  the  difference 
between  the  original  and  the  reconstructed  XBT  using  80  modes 
to  be  negligible,  of  the  order  of  10“*  degrees  Celsius, 
whereas  Figure  48,  which  shows  the  difference  between  the 
original  and  the  reconstructed  XBT  using  4  modes,  shows  a 
much  larger  overall  error  of  0.44  degrees  Celsius.  Table  III 
gives  the  mean  square  error  for  a  selected  number  of  modes. 
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C.  RECONSTRUCTION  AT  CAST  SITES 


In  addition,  the  objective  analysis  was  performed  at  cast 
sites  taken  within  the  10%  error  variance  contour  line  of 
Figure  32.  Thus,  instead  of  the  objective  analysis  being  done 
on  a  regular  grid,  the  procedure  reconstructed  synthetic 
XBTs  at  the  same  sites  where  the  original  XBTs  had  been 
taken.  This  was  done  as  a  check  to  ensure  that  the 
reconstructed  error  map  was  consistent  and  to  gain  a  measure 
of  how  much  error  there  was  between  the  original  XBTs  and  the 
synthetic  reconstruction.  A  selection  of  these  XBTs  are  shown 
in  Figures  50  -  54,  along  with  a  graph  of  their  associated 
RMS  error.  The  position  of  the  original  casts  can  be  found 
from  Figure  49. 

The  RMS  error  at  each  of  the  80  depth  setting  is  computed 
as  a  percentage  of  the  temperature  value  compared  to  the 
reconstruction  using  4  modes.  An  average  error,  expressed  as 
a  percentage,  is  then  obtained  for  each  XBT,  and  the  results 
are  averaged  over  the  set  of  XBTs  used  in  the  analysis.  The 
overall  error  between  the  synthetic  XBTs  compared  to  the 
reconstruction  using  the  four  original  EOFs  was  5%,  well 
within  the  10%  boundary. 

Reconstruction  of  all  15C  XBTs,  using  only  four  modes, 
gave  an  error,  when  compared  to  the  original  XBTs,  of  between 
6-7%.  The  overall  error  between  the  OA  reconstructions  and  the 
original  XBTs  was  found  to  be  between  10  and  11%. 
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D.  RECONSTRUCTION  AT  SELECTED  GRID  POINTS 

A  selection  of  XBTs  were  reconstructed  at  grid  point  sites 
and  the  error  compared  to  the  originals  that  were  likewise 
taken  at  the  same  points.  These  plots  and  the  associated 
error  graphs  are  shown  in  Figures  55  -  62.  The  position  of 
each  XBT  is  shown  in  Figure  63.  XBTs  71  and  88,  shown  in 
Figures  60  and  61,  are  displaced  from  the  nearest  grid  point 
(38.5  N  70  W) ,  to  which  they  are  compared.  These  profiles  are 
included  to  show  the  wide  range  of  variability  that  exists 
within  short  spacial  distances. 
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Figure  33  Reconstructed  XBTs  at  positions  37.5  N  72  W 
(a) and  37.5  N  71.5  W  (b) . 
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Figure  37  XBTs  766A  (a) 


Figure  39  XBTs  7144A  (a)  and  7143A  (b) 
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LATITUDE 


Figure  44  XBTs  7103A  (a)  and  761A  (b) 
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Figure  45  Positions  of  XBTs  shown  in  Figures  33-44 


Figure  46  XBT  7105A  showing  original  and  reconstructions 
using  4  inodes  and  80  inodes. 


Figure  47  Graph  showing  difference  between  original 
temperature  values  and  those  produced  through 
reconstruction  of  the  XBT  using  all  80  EOF's. 
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Figure  48  Difference  between  original  temperature  and  that 
produced  by  reconstruction  of  XBT  using  only  4  inodes. 
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Figure  49  Sites  used  for  selected  reconstructions- 


tetB^Celcha 


perceodfe  oiw 


Figure  50  a,  XBT  3  compared  with  reconstructions  using  4 
and  80  EOFs,  the  OA  being  performed  onto  the  site  of  the 
XBT  cast,  b,  rms  error  between  the  OA  and  the  4  and  80 


mode  reconstruction  for  all  depths. 
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Figure  51  a,  XBT  6  compared  with  reconstructions  using  4 
and  80  EOFs,  the  OA  being  performed  onto  the  site  of  the 
XBT  cast,  b,  rms  error  between  the  OA  and  the  4  and  80 
mode  reconstruction  for  all  depths. 


the  XBT  cast-  b,  rms  error  between  the  OA  and  the  4  and  80 


mode  reconstruction  for  all  depths. 
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Figure  53  a,  XBT  37  compared  with  reconstructions  using 
4  and  80  EOFs,  the  0A  being  performed  onto  the  site  of 
the  XBT  cast,  b,  rms  error  between  the  OA  and  the  4  and  80 
mode  reconstruction  for  all  depths. 


Figure  54  a,  XBT  63  compared  with  reconstructions  using 
4  and  80  EOFs ,  the  OA  being  performed  onto  the  site  of 
the  XBT  cast,  b,  rms  error  between  the  OA  and  the  4  and  80 
mode  reconstruction  for  all  depths. 


Figure  55  a,  XBT  14  compared  to  4  and  80  modes  and  to  OA 
reconstruction,  b,  RMS  error  between  OA,  4  modes  and  the 
original  for  each  depth. 
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Figure  56  a,  XBT  37  compared  to  4  and  80  modes  and  to  OA 
reconstruction;  b,  FMS  error  between  OA,  4  inodes  and  the 
original  for  each  depth. 


reconstruction;  b,  RMS  error  between  OA,  4  modes  and  the 


original  for  each  depth. 
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Figure  58  a,  XBT  121  compared  to  4  and  80  modes  and  to  OA 
reconstruction;  b,  RMS  error  between  OA,  4  modes  and  the 
original  for  each  depth. 
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Figure  59  a,  XBT  122  compared  to  4  and  80  modes  and  to  OA 
reconstruction;  b,  RMS  error  between  OA,  4  modes  and  the 


original  for  each  depth. 
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Figure  62  a,  XBT  151  compared  to  4  and  80  modes  and  to  OA 
reconstruction;  b,  RMS  error  between  OA,  4  modes  and  the 
original  for  each  depth. 
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Figure  63  position  of  selected  XBTs  casts.  Synthetic  XBTs 
were  calculated  at  grid  points  corresponding  to  original 
XBT  sites. 
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VI  THE  RESULTS 


A  succession  of  error  maps  using  reduced  data  are  shown  in 
Figures  64-70.  It  can  be  seen  that  the  data  were  taken  in  two 
natural  clusters,  in  the  discussion  that  follows,  only  the 
Western  cluster  is  considered.  Within  this  cluster,  most  of 
the  XBTs  lie  within  an  area  covered  by  the  10%  contour.  After 
reducing  the  total  number  of  XBTs  to  40,  the  area  covered  by 
the  30%  contour  is  still  on  the  order  of  the  size  of  the  area 
covered  by  the  original  10%  contour. 

The  experiment  was  refined  to  identify  a  specific  area 
that  lay  within  the  original  10%  error  contour  line. 
Additionally,  only  XBTs  taken  in  and  around  the  designated 
area  were  included  in  the  subsequent  analysis.  The  cluster  to 
the  east  was  removed,  plus  a  few  XBTs  laying  in  the  extreme 
north  of  the  analysis  area.  This  resulted  in  133  XBTs  being 
used  for  the  analysis  (Figure  71) .  The  aim  of  the  experiment 
was  to  reduce  the  data  set  until  the  30%  contour  intruded  into 
the  specified  area.  The  sequence  is  shown  in  Figures  72  - 
75.  The  30%  contour  crosses  the  borders  of  the  designated  area 
when  the  data  set  is  reduced  to  69  XBTs. 

It  was  concluded  that,  for  a  Gulf  Stream  meander,  a 
minimum  of  69  XBTs  is  required  to  adequately  reproduce 
synthetic  vertical  temperature  profiles  with  an  acceptable 
error  variance  of  30%. 
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Figure  64  Reconstruction  error  variance  using  100  out  of 
156  XBTs.  The  contour  interval  is  0.1 


Figure  65  Reconstruction  error  variance  using  50  out  of 
156  available  XBTs. 
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latitude  North 


Figure  68  Reconstruction  error  variance  using  25  out  of 
156  available  XBTs. 


Figure  69  Reconstruction  error  variance  using  20  out  of 
156  available  XBTs. 


Figure  70  Reconstruction  error  variance  using  10  out  of 
156  available  XBTs. 
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Figure  72 

Reconstruction  error  variance  using  100  out  of  the 
possible  133  XBTs. 


Figure  73  Reconstruction  error  variance  using  75  out  of 


the  133  available  XBTs. 
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Figure  74  Reconstruction  error  variance  using  70  out  of 
the  133  available  XBTs. 


Figure  75  Reconstruction  error  variance  using  68  out  of 
133  available  XBTs.  Note  the  30%  error  variance  line  which 
crosses  into  the  designated  area. 


VI  DISCUSSION 


A.  THE  RECONSTRUCTION 

The  difference  between  an  original  XBT  profile  and  its 
reconstruction  using  only  the  first  four  EOFs  has  been  of 
concern  throughout  this  study.  The  analysis  showed  the 
average  difference  was  of  the  order  6%.  Thus,  >efore  the 
objective  analysis  is  undertaken,  a  degree  of  error  has 
already  been  introduced  with  the  best  that  can  be  hoped  for 
being  a  contour  positioned  on  the  error  variance  map  accurate 
to  within  plus  or  minus  6  %.  As  a  result,  a  synthetically 
reproduced  XBT  will  have  associated  with  it  error  due  to  the 
objective  analysis  and  error  due  to  the  use  of  a  truncated 
series  of  four  EOFs.  However,  the  first  four  modes  account  for 
over  98%  of  the  variance,  and  reflect  a  minimum  number  that 
could  reasonably  be  used.  If  higher  accuracy  was  required, 
then  more  modes  could  have  been  considered  but  at  the  risk  of 
of  including  noise  from  individual  observations. 

All  these  sources  of  variability  are  included  in  an  80 
modes  solution.  The  current  situation  is, in  effect,  a  trade 
off;  loosing  some  of  the  fine  structure  due  to  the  small  XBT 
data  set  and  analysing  only  a  limited  number  of  vertical 
modes.  Nevertheless,  the  study  shows  that  a  limited  number  of 
vertical  modal  amplitudes  may  be  interpolated  using  objective 
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analysis  to  synthetically  create  XBTs  at  any  given  point 
within  the  region,  with  a  definitive  statement  as  to  the 
level  of  confidence  that  can  be  placed  in  the  reconstruction. 

B.  THE  NUMBER  OF  XBTs 

The  last  set  of  data  runs  in  this  study  (Figures  71-75) , 
provide  an  example  of  a  realistic  military  or  scientific 
scenario.  The  question  asked  was  how  many  XBTs  need  to  be 
taken  in  a  Gulf  Stream  meander  for  a  reasonable  estimate  of 
the  ocean's  vertical  temperature  structure  can  be  inferred  any 
where  within  the  meander? 

The  area  initially  chosen  was  the  area  within  the  10% 
error  variance  contour  and  reflects  the  area  that  was  most 
heavily  surveyed.  The  XBTs  surrounding  the  area  were  also 
included,  as  they  were  considered  to  represent  XBTs  that  would 
be  dropped  by  units,  whether  by  ship  or  aircraft,  that  were 
proceeding  to  or  away  from  the  area.  Overall,  XBTs  are  not 
dropped  at  regularly  spaced  intervals,  and,  although  not 
random  in  nature,  they  tend  to  reflect  a  distribution  that 
would  be  expected  to  be  produced  by  several  surface  units 
attempting  to  track  a  covert  submarine. 

The  area  noted  in  Figures  71-75  is  approximately  1400 
square  miles  and  was  initially  surveyed  by  133  XBTs.  The 
analysis  indicates  that,  given  a  confidence  level  of  30% 
error,  the  same  area  could  have  been  adequately  sampled  by  70 
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XBTs.  This  is  a  saving  of  nearly  50  percent  in  XBTs,  but,  more 
importantly,  this  study  indicates  that  an  effective  analysis 
can  be  achieved  in  a  complicated  region  with  relatively  few 
XBTs.  Although  this  study  has  utilized  data  from  within  only 
one  Gulf  Stream  meander,  it  provides  a  general  indication  of 
the  amount  of  observations  that  would  be  needed  within  other 
Gulf  Stream  eddies  or  meanders. 


C.  OPTIMAL  SPACING 

The  determination  of  the  spacial  correlation  matrices 
resulted  in  parameter  b,  the  e  folding  distance,  to  be  defined 
and  calculated  for  each  of  the  modal  amplitudes.  This  distance 
places  a  limit  on  the  separation  that  can  exist  between  two 
observations  to  be  included  in  the  analysis.  From  Table  II  it 
can  be  seen  that  the  second  modal  amplitude  gives  the  smallest 
value,  a  distance  of  25  km.  This  value  represents  the  maximum 
distance  of  separation  that  should  exist  between  two  adjacent 
observations.  For  the  purpose  of  economy  and  military 
logistics  the  figure  represents  the  optimal  spacing  that 
should  exist  between  XBT  cast  sites. 

The  value  of  25  km  is  approximately  half  that  of  the 
Rossby  radius  of  deformation  and  is  suggestive  that  a  smaller 
grid  scale  would  have  been  more  appropriate. 


97 


0.  RECOMMENDATIONS 


To  valididate  the  claim  that  25  km  is  a  good  optimal 
distance,  it  would  be  of  value  to  extend  the  study  to  consider 
regularly  spaced  XBTs  (generating  them  synthetically,  as  the 
data  from  any  real  survey,  by  its  very  nature,  will  tend  to 
have  been  erratically  sampled  data  in  terms  of  both  time  and 
space)  with  the  distance  between  adjacent  casts  being 
gradually  extended  until  the  resulting  error  variance  becomes 
unacceptable. 

The  current  analysis  also  does  not  take  into  account  the 
fact  that  each  XBT  was  taken  at  different  times.  It  was 
assumed  throughout  the  study  that  all  XBTs  were  valid  at  the 
analysis  time.  The  study  could  be  extended  to  take  time  into 
account,  with  the  interpolation  being  adjusted  to  allow  for 
an  optimal  value  to  be  chosen  both  in  terms  of  time  and  space 
(see  Carter  1982). 

A  different  correlation  function  could  also  have  been 
fitted  to  the  cross  flow  and  along  flow  directions.  This  has 
value  as  it  helps  to  account  for  the  rapid  changes  that  take 
place  across  the  Gulf  Stream  front  as  opposed  to  the  expected 
similarity  in  values  taken  along  the  front.  In  this  study,  the 
casts  were  taken  within  a  well  developed  horseshoe  shaped 
meander  so  it  was  decided  to  assume  homogeneous  statistics 
using  the  same  correlation  function  in  all  directions. 
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However,  the  use  of  a  non  isotropic  field  should  be 
considered. 

A  major  extension  to  the  study  would  be  to  obtain  the 
principle  modes  by  including  data  from  other  Gulf  Stream 
eddies  and  meanders  so  as  to  build  a  climatology  of  Gulf 
Stream  eddies.  It  is  likely  that  the  modal  decomposition  of 
a  projection  matrix  defined  from  a  larger  data  set  would 
remove  the  spurious  effects  evident  in  the  current  study  and 
allow  for  an  improved  reconstruction  of  the  data  when  using 
the  four  principle  modes.  It  is  considered  that  a  climatology 
of  eddies  rather  than  a  climatology  of  the  North  West  Atlantic 
would  be  of  greater  value  in  attempting  to  empirically  model 
XBTs  within  the  Gulf  Stream  region. 

It  is  noted  that  the  surface  layer  is  poorly  modelled, 
suggesting  that  two  analyses  may  be  required.  One  analysis  for 
the  surface  layer,  the  upper  80  metres,  and  the  second  for 
the  deeper  water  below  the  thermocline. 
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VII  CONCLUSION 


The  creation  of  synthetic  XBTs  at  regular  locations  within 
a  Gulf  Stream  meander  by  the  use  of  an  objective  analysis  of 
modal  amplitudes  produced  from  the  decomposition  of  the 
vertical  temperature  correlation  projection  matrix  has  been 
shown  to  be  of  value.  Although  there  is  a  degree  of  error  in 
the  reconstruction,  the  value  of  the  error  is  explicitly 
stated. 

Using  the  error  variance  field  generated  from  the 
objective  analysis,  it  has  been  shown  that  within  a  1400 
square  mile  region  of  a  warm  Gulf  Stream  meander  a  minimum  of 
69  XBTs  need  to  be  taken  in  order  for  a  synthetically  produced 
XBT  to  be  within  30%  of  its  true  value. 

The  spacial  correlation  statistics  indicate  that  the 
optimal  distance  between  XBT  cast  site  must  be  25  km  or  less. 
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Oct  2 
Oct  4 
Oct  4 
Oct  4 
Oct  4 
Oct  4 
Oct  4 
Oct  4 
Oct  4 
Oct  5 
Oct  5 
Oct  5 
Oct  6 
Oct  6 
Oct  6 
Oct  6 
Oct  7 
Oct  7 
Oct  7 
Oct  7 
Oct  7 
Oct  9 
Oct  9 
Oct  9 
Oct  9 
Oct  9 
Oct  9 
Oct  9 
Oct  9 
Oct  10 
Oct  12 
Oct  12 
Oct  12 
Oct  12 
Oct  12 
Oct  12 
Oct  12 


88275:03:03 

?,7 

9.1 

N 

71 

50.5 

W 

88275:03:57 

37 

10.0 

N 

71 

59.9 

W 

88275:07:36 

37 

38.8 

N 

71 

44.8 

w 

88275:08:57 

37 

54.1 

N 

71 

31.1 

w 

88275:10:23 

37 

58.4 

N 

71 

10.2 

w 

88275:11:27 

38 

0.8 

N 

70 

51.0 

w 

88275:12:20 

37 

56.7 

N 

70 

45.9 

w 

88275:17:41 

38 

5.9 

N 

70 

29.8 

w 

88275:20:18 

38 

28.0 

N 

70 

29.9 

w 

88275:21:14 

38 

28.1 

N 

70 

17.6 

w 

88276:01:09 

38 

6.3 

N 

70 

15.0 

w 

88276:02:13 

38 

6.5 

N 

69 

59.3 

w 

88276:05:15 

38 

28.4 

N 

69 

59.8 

w 

88276:06:09 

38 

27.9 

N 

69 

45.2 

w 

88276:12:44 

38 

27.7 

N 

69 

30.0 

w 

88276:13:54 

38 

27.9 

N 

69 

14.3 

w 

88278:13:43 

39 

0.1 

N 

70 

43.5 

w 

88278:15:06 

39 

0.0 

N 

70 

30.0 

w 

88278:16:20 

39 

0.0 

N 

70 

14.9 

w 

88278:17:31 

39 

0.1 

N 

70 

0.5 

w 

88278:18:52 

38 

59.9 

N 

69 

43.6 

w 

88278:19:55 

39 

0.0 

N 

69 

30.1 

w 

88278:21:06 

39 

0.0 

N 

69 

15.0 

w 

88278:22:28 

39 

0.0 

N 

68 

58.1 

w 

88279:11:21 

88279:15:25 

88280:13:12 

88280:19:23 

88280:21:07 

88280:22:36 

88281:05:37 

88281:07:34 

88281:08:42 

88281:11:06 

88281:11:42 

88283:11:44 

88283:13:43 

88283:14:26 

88283:15:25 

88283:15:35 

88283:16:40 

88283:17:40 

88283:19:00 

88284:15:39 

88286:11:11 

88286:12:35 

88266:13:51 

88286:15:11 

88286:16:28 

88286:19:35 

88286:20:46 


45.4 

44.7 
45.0 
30.0 

39.6 

24.1 
29.9 

10.8 

18.1 

22.1 

18.4 

20.4 

30.3 
11.1 
4.5 
11.8 
11.8 
2.4 
0.0 

50.0 

17.6 
30.0 
30.0 
30.1 
30.0 
29.8 

44.7 

58.4 


xbi  li; 
xbt  li; 
xbt  1 1< 
xbt  in 
xbt  lr 
xbt  HI 
xbt  1H 
xbt  12C 
xbt  121 
xbt  12  j 
xbt  123 
xbt  124 
xbt  125 
xbt  126 
xbt  127 
xbt  129 
xbt  130 
xbt  131 
xbt  132 
xbt  133 
xbt  134 
xbt  135 
xbt  137 
xbt  138 
xbt  139 
xbt  140 


70 

30.0 

W 

xbt 

141 

71 

9.3 

W 

xbt 

142 

71 

9.3 

W 

xbt 

143 

71 

8.6 

W 

xbt 

144 

70 

19.3 

W 

xbt 

145 

70 

28.1 

W 

xbt 

146 

70 

12.3 

W 

xbt 

147 

70 

0.3 

W 

xbt 

148 

69 

50.6 

W 

xbt 

149 

69 

29.3 

W 

xbt  150 

69 

17.1 

W 

xbt 

151 

69 

9.4 

W 

xbt 

152 

68 

58.4 

W 

xbt 

153 

68 

56.1 

w 

xbt 

153 

68 

45.4 

w 

xbt 

154 

68 

30.0 

w 

xbt 

155 

68 

15.0 

w 

xbt 

156 

69 

23.5 

w 

xbt 

1 57 

68 

0.1 

w 

xbt 

158 

68 

20.6 

w 

xbt 

159 

68 

40.0 

w 

xbt 

160 

69 

0.3 

w 

xbt 

161 

69 

20.2 

w 

xbt  162 

69 

30.0 

w 

xbt 

163 

69 

29.9 

w 

xbt 

164 

c 


APPENDIX  IB 


program  loadbathys 
c 
c 

c444444444444444444444444444444444444*4444444444#44444444444444444444444 

c  this  file  loads  the  bathys  into  an  array,  and  calculates  correlation. 

£.4444*4444  +  44*4444444444444444444444444444444*444444*4444444444444444444 


integer  m,»,p, iy, iz,q 

real  sumab,  corro,  a,  b,  y ,  z  ,  tempvar  (80,  156) .  sumvar  (80) 

integer  counta(80, 156) , count 

real  deptha(30, 156) .corrol (80,82) 

real  tempa (80,156), mean ( 82 ) , times, volt , qual 

rea  1  sumsqa ,  eumsqb,  dept he  ( 2000 ,  156 )  ,  tempo  (2000,  156) 
character  latMO,  long*12,  t  ime*4 
character  recnuma (80, 156) *8, remim(156) *9 
integer  yearday 


c  t************************************************************** 


c  LOAD  IN  DATA 


c 


c 


210 

220 

c 


write(*, * 
open (unit 
open (unit 
open (unit 

p  r  0 


'loading  bathys' 

4,  file  =  'deepname ', status 

20,  f i]«  =  'name' ) 

21,  file  =  ' namepos ' ) 


'old' ) 


do  200  n=l, 156 

read (4 , ' (a9 ) ' , end 
P  =  P  +  1 
write(20, *)  p, ' 
open(unit=  3, file 


=  2  30)  renum(ti) 

’ , renum(n) 

= ' /usr/whitney_dl/xbt/ ’ //renum(n) ) 


rewind  3 
read( 3,  *) 

read (3, 210)  yearday, time 
read (3, 220)  lat. ,  long 
read ( 3 , * ) 
read ( 3 , * ) 

wr i t e ( * , * )  yea r day .time 
write(*,*)  lat, long 
write (21,*)  lat, long 
format (9x, i3,t29,a6) 
format (6x,al0,t25,al2) 


do  240  m  =  1,20.00 

read (3,*, end  =  200)  t imes , depthc (m, n ) , 


$ 


volt, tempo (m, n) , qual 


240 

cont inue 

write ( * , * 

)  P 

200 

continue 

230 

close(4 ) 
close ( 3 ) 

c 

write ( * , *)  p, ' 

bathys  loaded' 

**»**************»*************»********»**»»*«*********«»***** 

c  This  section  calls  redata  and  calculates  temp  at 
c  10m  increment  depths,  starting  at  5m,  Cor  each  bathy  and  loads 
c  them  into  arrays.  Also  finds  mean  temp  for  given  depth  taken  over 
c  all  bathys. 

— *«»***********«********«***»**«**«*4*****»«**«*«***«**«*»****4** 

c  CALCULATE  SMOOTHED  TEMP  AMD  MEAN  FOR  GIVEN  DEPTH 
c 

open(unit  =12  .file  =  'meantemp' ) 

do  40  y  =  5,800, 10 
sum  =  0 

iy  =  1* ( (y-5)/10) 

call  redata (y, iy, remim, count , count#, recnuma, depths, 
$  tempa, depthc, tempc, p) 

c 

do  30  n  =  l,p 

sum  =  sum  +  tempa (iy,n) 

30  continue 

mean(iy)  =  sum/count, 
write! 12,*)  y.menn(iy) 

40  continue 


c  SEND  TO  OUTPUT 
c 

open(unit  =  13  ,  file  =  'prof i le .mat ' ) 

do  100  n  =  l.p 

do  110  Iy  =  1,80 
writ® (13,*)  tempa (iy, n ) 

110  continue 

100  continue 

c  now  loop  through  each  depth  calculating  the  correlation  compared 
c  with  the  shallow  depth. 

e************************************************ •*******«***•«• 

c 

c  CALCULATE  CORRALATION 

ccc  open (unit  =  10, file  =  'lbokat') 

open(unit  =  8, file  =  'corrl.mat’) 


•  the  "shallow*  depth  <  the  m  loop) 


c 


c 


500 

cccc 

cccc 


510 


c 

310 


305 

c 

300 

c 


do  300  y  =  5, 800, 10 

iy  =  1  +  <y-5) / 10 
write!*,*)  'depth  y  , iy 

do  305  z  =  5,800,10  !  the  "deep*  depth  (the  n  loop) 

!  compare  each  deep  with  shallow 
set  constants,  counter  etc  to  zero 
a  =  0 
b  =  0 
sumab  =  0 
sumsqa  =  0 
sumsqb  =  0 
corro  =  0 
iz  =  1+ ( z-5 ) / 10 

do  310  n  =l,p  !  loop  through  each  "deep* 

m  =  n  !  record 

q  =  n 

if  { (d«ptha ( iy , m) ) . eq.  (0.0) )  goto  310 
if ( (dept ha < iz,  n) ) . eq. ( 0 . 0 ) )  goto  310 

if  ( r ecnumn ( iy , m) . eg. r ®cnuma ( i z, n ) )  then 

m,  iy , dept  ha ( iy , n) , tempa ( iy , m) , recnuma ( iy , m) 

n,  iz, dept ha ( iz, n) , tempa { iz,m) , recnuma ( iz, n) 
a  =  tempa (iy,m)  -  mean(iy) 
sumsqa  =  sumsqa  +  a* *2 
b  =  tempa (iz.n)  -  mean(iz) 
sumsqb  =  sumsqb  +  b**2 
sumab  =  sumab  4 (n*b) 

else 

do  510  m  - 1 , p 

if  (recnuma  ( iy,  ni)  .  eq.  recnuma  ( iz,  n) )  goto  500 
font inue 
goto  310 

end  i  f 

continue  !  with  next  record 

corro  =  sumab/sqrt (sumsqa  ‘sumsqb) 
corrol ( iy, iz)  -  corro 
write(8,‘)  corrol ( iy, iz) 

!  next  “deep"  depth 

!  nekt  "shallow*  depth 


continue 

continue 

end 


write( 10, * ) 
write(l0, *) 


c 


APPENDIX  2B 


subroutine  redata (z, iy, renum, count, counta, recnuma, deptha, 

$  tempa, depthc, tempc, p) 

c  this  program  finds  temp  for  a  given  depth  for  all  records 
c  and  stores  the  value  of  temp  and  depth  in  arrays  and  passes  them 
C  back  to  vertcorro. f  .  The  depth  for  given  temp  is  found  by 
c  averaging  over  a  10m  bin.  The  value  of  depth  used  is  passed  in 
c  the  call  from  main  program. 

£»*************»******+*+*»*********»*********»****»************** 


parameter  (mm  -2) 

real  t, depth (mm) , temp (mm) , z,  tempo (2000,  156) 
real  deptha (80, 156) , tempa (80, 156) , depthc(2000, 156) 
character  renum ( 156) *8 , recnuma (80 , 156 ) * 8 
integer  n, count , counta (80, 156) , iy 
integer  counter , nodatapt , p 
real  add 

{.ntmttutttnttmiiHiiHiiutiiHiittttiitiHKiiuuHiHit 

c  open (unit  =  10, file  =  'datapt')  !  testing  for  data 

c 

e  SET  COUNT  the  number  of  records  processed  for  a  given  depth 
c 

count  =  0 
c 

1010  do  600  ns  l,p  !  loop  through  each  bathy 


c 

c 


c 


c 

c 

c 

c 


add  =  0  !  the  sum  of  data  points  ut  d  in  a  bin 

counter  =0  !a  count  of  number  of  data  points 
!  used  in  a  bin 

do  2000  m  s  1,2000  !  loop  through  all  data  points 

dept’i(l)  =  depths  (m,  n) 
temp(l)  -  tempc(m,n) 

SELECT  DATA  POINTS  IN  BIN 

if  (depth!) ) .ge. z-5.and.depth( 1) . le. Zf5)  then 

counter  =  counter  +1  (  add  up  number  of 

data  points 

add  =  add  +  temp(l)  '■  sum  values  of 

!  data  points 

t  a  add/counter  !  mean  temp  for 

!  given  depth  for 
!  given  bathy 


if  (depthc (nu  1 , n) . gt . 74 5)  goto  5000 
if  ( ( t.empc (m4  1, n) )  .  le.  (0 . 0) )  goto  5000 
goto  2000 
else 

goto  2000 
C 

c  LOAD  VALUES  INTO  ARRAYS  FOR  PASSING  BACK. 


c 


-500"  count  =  count  +  1  !  increment  counter 

count* (i  n)  =  count 

t*mpa(iy,n)  =  t  *  value  of  t  placed  in  array 
recnuma ( iy, n) =  renum(n)  !  name  of  record  beng  read 
depth* ( iy, n)  =■  z 

nodatapt  =  counter *  count  !  calculate  number  of 
c  !  data  points 

goto  600 
end  if 

2000  continue 


600  continue 

c  write { 10, *) z, nodatapt  , count , counter  !  sends  to  file 

END 


c  APPENDIX  3B 

c 

program  distcorro 


parameter (max  *  56) 
implicit  real*8  (a-z) 
real* 8  dpr  /57 . 29577951308232/ 
real*4  dist (max, max) 

real  declat, declong, lat (max) , long(max) .maxvalue, noofbins 
integer  n,m, z, lati , longi, p, width 
real  mode(4,max)  ,corr(4,max)  ,  value  (max.  max) 
integer  inoofbins 

c  ************************************************************ 

c  This  part  calculates  the  distance  between  any  two  points  and  stores  in 
c  an  array 


open (unit  =  3, file  =  'decpos .mat ' ) 
open (unit  =  1,  file  =  'deeppos') 


do  10  n  =  l.max 

read(l, 20)  lati,  declat , longi , declong 
20  format (x,i2,x,f5.2.4x,i3,x,f5.2) 

declat  =declat/60 
declong  =  declong  /SO 
lat (n)  =  lati  +  declat 
long(n)=  longi  -  declong 
write (3,*)  long (n) , lat (n) 

10  continue 

c 

do  30  n  a  l.max 
do  40  m  = l.max 


c 

t  =  (long(n)  -  long (m) ) /dpr 
dsinr=  dsin(t)  *  dcos ( lat (m) /dpr ) 
r  =  dasin(dsint) 

dsins  =  dsin( lat (m) /dpr ) /dcos (r) 
s  =  das  in (dsins) 

dcosrl  =  dcos  (r)  *  dcos(s  -  (lat(n)  /dpr)) 
d  =  dacos(dcosd)  *  dpr 
c 

dist(n.m)  =  d  *  111.6114 


c 

40  continue 

30  continue 


c 

c 


c  Now  load  the  modal  amplitudes,  (just  the  principle  ones  for  now) 


open(unit  =  2,  file  =  'single. mat' ) 


do  110  z  =  1,4 
do  100  n  =  l,max 

read(2,*)  mode(z.n) 

100  continue 

110  continue 

c 

c  ***************************************************** 

c 

c  find  max  distance  for  this  data  set 
c 

tnaxvalue  =  0.0 
do  210  n  =l,max 
do  200  m  =  l,max 

if(  dist (n,m) .gt .  maxvalue)  maxvalue  =  dist(n.m) 
200  continue 

210  continue 

write(*,M  'maxvalue  =  maxvalue 


c  create  bins  by  distance  . 
c 

write(*,*)  'enter  bin  width' 
read(*,*J  width 
end  =  max  *  max 
noofblns  =  maxvalue /width 
inoofbins  =  int (noofbins) 


c 

c  Calculate  correlation  as  function  of  distance, for  each  bin(p). 


open)  unit.  =  10,  file  ='smith.l') 
open?  unit  =  11, file  ='smith.2') 
open(  unit  =  12, file  = 'smith. 3') 
open(  unit  =  13, file  =  'smit.h.4’) 

do  700  z  =1,4 
do  600  p  =  1, inoofbins 

call  calc(mode, z , p, cor r , width, dist ) 

continue 

********»•**»»»*»«*»»*»»•» 

send  correlation  for  each  bin  to  file  smith 

do  400  p  =  1, inoofbins 
write<9  ♦  z, M  p*width,corr iz.p) 


600 

c 

c 

c 

c 

c 


400 

700 

c 


continue 

continue 


i  *****  *  i 


h«+4********+**+**«***«*****#**+**4*ft 


C 

©nd 


c 

subroutine  calc (mode, z, p, corr, width, dist) 
c 


parameter (max  -  56) 

real  modea  ,  modeb,  meantnodea,  cot  L  (4,  max) 

real  meanmodeb, summodea, summodeb, sumab, sumsqa, sumsqb 

real  mode ( 4 , max ), dist (max, max) 

integer  p, counter, n,m,a,b,z,width,ct  number b ( 5000 ) 
integer  numbers (5000 ) 


r 

c 

c 

c 

r* 

C 


open!  unit 
open (unit 
open (unit 
open (unit 


1, 

file  = 

' mean  * ) 

2, 

file  - 

’  look  *  ) 

3, 

file  - 

’  a  loop  * ) 

4, 

file  = 

’  bloop ' ) 

c 

c 


set  counter  and  all  variables 


back  to  zero  for  a  new  p  (bin) . 


counter  -  0 
modea  =  0 
modeb  -  0 
meanmod»n  r  0 
meanmodeb  -  0 
summodea  =  0 
summodeb  =  0 

SHIM  b  =  0 

sumsqa  0 
sumsqb  s  0 
corr(z,p)  =  0 


c 

c  Determine  which  data  points  are  used  for  a  given  bin. 


c 


HO 

300 


do  300  u  - 1 , ma x 
do  310  m  =  n, max 
a  =  p'width 
b  -  p*widt)i  -25 
wr  i  t »  <  * , * )  n , m 

if  (dist  (n,  m)  .  It .  a  .  and  .  dist  (n,  m)  .  cie  .  b)  then 


counter  -  counter  *1 
c  -  counter 
write)*,*)  a,c 
numbers  (c)  =■  n 
numberb(c)  =  m 


end  i  f 

cont inue 
continue 


c 


write(*,*)  p, p* width, c 


c  the  value  c  i.9  the  number  of  data  point  pairs  in  a  bin. 
c  numbera  is  the  array  number  containing  the  first  point  and 
c  numberb  is  the  array  number  containing  the  second  point . 
c 

c  now  add  up  value  of  all  data  points  used  for  given  bin. and 
c  calculate  the  mean  for  "a*  loop. 

counter  =  0 

do  400  n  =  l.max 
do  410  in  -  1 ,  c 

i f (n . eg . numbera (m) )  then 
summodea  =  summodea  +  mode(z,n) 
write (3,*)  p, n, m, mode ( z , n) 
counter  =  counter  +  1 
goto  400 

else 

goto  410 
end  i  f 
cont inue 
cont inue 

meanmodea  =  summodea /counter 

r.****«»**********************»****»*****»«*********»*4****** 

c  now  do  same  for  *b"  loop. 


410 

400 


counter  «  0 
do  500  n  =  l,max 
do  5 1 0  in  -  1 ,  c 

l  f  (n .  eg .  numberb(m) )  tlien 
summodeb  -  summodeb  +  mode(z.n) 
counter  =  counter  +  1 

c  writeM,*)  p,  n,  m,  mode  (m,  n ) 

goto  500 

else 

goto  510 
end  i  f 

510  continue 

500  continue 

meanmodeb  =  sununodeb/counter 


c  write (1,*)  z, p, reanmodea, meanmodeb 

r**********#********************l*****»****************4****** 

c 

r  Calculate  component!?  for  corrolation 

c 

do  320  n  =  l,  max 
do  330  m  = 

a  =  p* width 


b  =  p*width  -25 

if  (dist (n, m) . It . a . and . dist (n,m) . ge . b)  then 


modea  =  mode(z,n)  -  meanmodea 
modeb  =  mode(z,m)  -  meanmodeb 
surnab  =  surnab  *  modea *modeb 
sumsqa  =  sumsqa  +  modea* *2 
sumsqb  =  sumsqb  +  modeb* *2 


end  if 


330 

cont inue 

320 

continue 

c 

c  calculate 

correlation, 

for  passing  back. 

corr ( z,  p) 

=  (sumab) /sqrt (sumsqa*sumsqb) 

c 

write( * , 

*)  z,p, corr (z,  p) 

end 


c 


APPENDIX  4B 


c 

c  This  will  be  a  program  that  determines  correlation 
c  for  any  distance  by  fitting  to  the  data  in  'smith', 
c 

******************* 

parameter!!  =  8) 
integer  m, j , n, count , istep 

real  cr (4 , i) , cn (4 , 28 ) , error ( 4, i ) ,  suraerror , d(4 , 28) , olderror 

real  a, b, newerror , deltaa, deltab,deltaeal,deltaea2,r 

real  deltaebl , deltaeb2 , incrementa, incrementb 

real  suma, sumb, deda , dedb 

real  alfa (4) , beta ( 4 ) 

c  ***************♦***********»**+************»*****♦*****♦«*»* 

open(  unit  =  11,  file  =  'smith.!') 

open  (  unit  =  12,  file  •=  'smith. 2') 

open(  tinit  =  11,  file  =  'smith. 3') 

open (  unit  =  14,  file  =  'smith. 4') 

open(  unit  =  2,  file  *  'look') 

c 


sumerror  =  0 

c  ****»»*****♦*****♦***♦»»*»»»*»** 

c  load  in  values  from  ' smith. a  1 1 ' . 
c  give  initial  values  of  a  and  b 

alfa ( 1 )  =  130 
beta(l)  =  60 
alfa (2)  =  60 
beta ( 2 )  =  30 
alf a ( 3 )  =  45 
beta ( 3 )  =  20 
al  fa (4 )  =  50 
beta(4)  =  20 
c 

do  200  m  =  1,4 


c  reset  all  variables  to  zero 

olderror  =  0 
newerror  =  0 
deltaeal  =  0 
deltaea2  =  0 
deltaebl  =  0 
deltaeb2  =  0 
incrementa  =  0 
incrementb  =  0 
suma  =  0 


sumb  =  0 
deda  =  0 
d«db  b  0 
count  a  0 

c  load  in  the  modes 

write(2,*)  'i  the  number  of  bins  ='  ,  i, >******** 
write(2, *} 
do  10  n  =  1,22 

read <10  f  m,*)  d(m,n),  cn(m,n) 
c  write{*,*)  n,  cn(m,n) 

10  continue 

write(2, *) 


a  =  alfa(m) 
b  =  beta(m) 
r 

c 

c  Calculate  cr 

do  20  n  =  1 ,  i 

r  =  c1(m,n) 

cr (m, n)  =  (1  -  <v/a ) * *2) *exp ( - (r /b) * *2 ) 

20  continue 

c 

r  it****************************************************** 

c 

c  calculate  value  of  error  {  that  is  to  be  minimized) 
c 

do  30  n  =  1 , i 

error(m,n)  =  (cr(m,n)  -  cn(m,n))**2 
sumerror  =  sumer ror  +  error (m,n) 

30  continue 

olderror  =  sumerror/n 
c 

c 

c  ******************************************************* 

c 

c  calculate  the  delta  error, and  decide  whether 
c  to  add  or  subtract  the  increment 

deltna  =  0.001 
deltab  =  0.001 
c 

suma  =  0 
sumb  s  0 


do  50  n  =  1, i 

deda  =  (4* (r**2) *exp<- (r/b) **4) / (a* *3 ) ) * { 1- (r/a) * *2 ) 

$  -  4* (r**2) * (exp (-(r/b)**2))*cn(m,n)/(a**3) 

suma  =  suma  +  deda 
50  continue 

deda  =  suma/n 

deltaeal  =  deda  *  deltaa 

deltaea2  =  deda  *  (-  deltaa) 

if  (deltaeal . It .deltaea2)  then 
incrementa  =  deltaa 
else 

incrementa  =  -deltaa 
end  if 

c  write(*,*)  deda, deltaeal , deltaea2 , incrementa 

do  60  n  =1 , i 

dedb  =  (4*(r*M)Mexp(-(r/b)**4))/b**5) 

$  *  (l-2MJ./a)**2  +  (r/a )  *  *  4  ) 

$  +  (4*(r *  *2) *cn (m, n) *  (exp( -  (r/b)**2))/b**3) 

$  *  ( ( r /a ) *  *  2  -1) 

sumb  =  sumb  +  dedb 
60  continue 

dedb  =  sumb/n 

deltaebl  =  dedb  *  deltab 

deltaeb2  =  dedb  *(-  deltab) 

if  (deltaebl . It .delta°b2)  then 
incrementb  =  deltab 
else 

incrementb  =  -deltab 
end  if 

c  write{*,*)  dedb, deltaebl, deltaeb2, incrementb 

c 

c 

c  ♦♦♦**+**♦*******#********♦*♦****«#**♦****♦*■»****«*♦♦**»***»♦♦ 
it*********************************************************** 

c 

c  increment  a  and  b. 

c 

100  count  s  count  f  1 

a  =  a  +  incrementa  !  normally  ♦  -  -  + 

b  a  b  -  incrementb  !  normally  -  -  +  + 

c 

r  *********************************************************** 

c 


c  calculate  cr  again 
c 

do  70  n  =  1, i 
r  =  d(m,n) 

cr(m,n)  =  (1  -  <r/a) **2) *exp(- (r/b) **2) 

70  continue 

c 

Q  **************************************************** 

c 

c  calculate  error  again,  with  new  a  and  b 
c 

sumerror  =  0 
do  80  n  =  1, i 

error(m,n)  =  (cr(m,n)  -  cn(m,n))**2 
sumerror  =  sumerror  +  error (m,n) 

80  continue 

newerror  =  sumerror /n 

write(*, ♦)  m, count, newerror , a,  b 
c  stop 

£  *****************************1  *********************** 

c 

iE  (newerror .gt .olderror)  then 
write(2, * ) 

write(2,*)  'for  minimized  error,  mode'.m 
write(2,*)'  a  b  error 

write(2,*)  a, b, newerror , count 

write(2, *) 

goto  195 
else 

olderror  =  newerror 
goto  100 
end  i  f 

195  alfa(m)  =  a 

beta(m)  =  b 
200  continue 

call  corro(alfa, beta, i) 
close (11) 
close (12) 
close(13) 
close(14) 
end 


iterations' 


subroutine  corro(alfa, beta, i) 
parameter (i  =  8) 
integer  m,n 

real  alfa(4) , beta (4) , a, b, cr (4 , i ) , bin, r 


open (unit  =  11, file  ='corl') 
open (unit  =  12, file  ='cor2') 
open(unit  =  13, file  ='cor3') 
open(unit  =  14, file  ='cor4') 


c 

alfa(l) 

= 

115.37 

c 

beta (1) 

a 

90.37 

c 

al fa (2) 

a 

84.18 

c 

beta (2) 

a 

59.18 

c 

alfa(3) 

= 

145.46 

c 

beta ( 3 ) 

a 

45.42 

c 

alfa(4) 

a 

51.86 

c 

beta (4) 

a 

46.86 

bin  =  25 
c 

o  calculate  cr 
c 

do  60  m  =1,1 
a  =  alfa(m) 
b  =  beta(m) 

write(*,*)  alfa (m) , beta (m) 

do  7  0  n  =  1 ,  i 
r  =  n  *  bin 

cr(m,n)  a  ( 1  -  ( r /a ) **2) *exp (- (r /b) *  *  2) 

c  print  out 

write(10  +  m  ,*)  r,cr(m,n) 

70  continue 

60  continue 


end 


c  APPENDIX  SB 


C  SPACE-TIME  OBJECTIVE  ANALYSIS/STATISTICAL  FORECAST  PACKAGE 
C  USING  GENERALIZED  OBJECTIVE  ANALYSIS  ROUTINE 
C 

C  (C)  COPYRIGHT  EVERETT  CARTER  1981 

C 

C  uses  NCAR  double  precision  matrix  inverter  INVMTX 

C 

C  UPDATES! 

C 
C 
C 
C 
C 
C 
C 

C  IER  is 

C  =0 

C  >0 

C  =-l 

C  =  -3 

C 

c 
c 
c 

c  THE  MAIN  PROGRAM  MUST  SET  UP  THE  DIMENSIONS  AS  FOLLOWS 

C  (FOR  A  33X33  FIELD) 

C  DIMENSION  PSK1089)  , XOBS ( 1089 , 2) , TOBS ( 1089 ) 

C  DIMENSION  X ( 2 , 1089 ) , THETA (1089) , EPS ( 1089 ) 

C  DIMENSION  PARS I (20) , T(20) , rARX(2,  20) 

C  COMMON  BLOCK  ERR  CONTAINS  THE  OBSERVATION  ERROR  PARAMETERS 
C  E  IS  THE  MEAN  SQUARE  NOISE  LEVEL  IN  TERMS  OF  PERCENT  OF  VAR 
C  COMMON / ERR / E 

C  THE  FUNCTION  F  IS  THE  CORRELATION  FUNCTION 
C  EXTERNAL  F 

C  M  IS  THE  TOTAL  NUMBER  OF  GRID  POINTS 
C  LIMIT  IS  THE  MAXIMUM  NUMBER  OF  INFLUENTIAL  POINTS 
C  DATA  LIMIT/ 10/ 

C  DATA  M/1089/ 

C  DATA  DIST.T1M/100. , 20. / 

C  PS I  THE  OBSERVATION  VALUES 

C  XOBS  THE  OBSERVATION  POSITIONS 

C  TOBS  THE  OBSERVATION  TIMES 

C  TCEN  THE  CENTRAL  INTERPOLATION  TIME  (PREDICTION  TIME) 

C  X  THE  INTERPOLATION  POSITIONS 

C  THETA  THE  INTERPOLATION  VALUE  OF  THE  COMPLETED  FIELD 

C  EPS  THE  INTERPOLATION  ERROR  LIMIT  OF  THE  COMPLETED  FIELD 

C  E=0 . 05 

C  EXAMPLE  MAIN  LOOC 

C  DO  150  KX= 1 , M 

C  CALL  SELECT(LIMIT,X(KX,  1.)  ,X(KX, 2)  .TCEN. XOBS, TOBS.  PSI, 

C  1  PARSI. PARX.T.N, NOBS, DIST, TIM) 


21  Aug  1984  —  Modified  package  so  that  GETAVE  is  called 

within  OBJ AN.  also  added  COMMON  block  CBLOCK 
iti  order  to  reduce  correlation  function  calls 
8  Aug  1984  —  added  routine  GETAVE,  to  remove  weighted  mean 
27  DEC  1983  --  expanded  IER  flags 
3  NOV  1983  --  added  poor  matrix  inversion  Warning 

an  error  flag  for  ob.tan 
for  no  eriors  detected 

for  matrix  inversion  errors  (see  matrix  inversion  routine) 
for  no  input  data  (a  Warning--  not  fatal) 

for  poor  matrix  inversion,  it  did  it  but  the  inversion  was 
nearly  NUMERICALLY  singular 


C  CALL  OBJAN ( TARSI , PARX, T, NOBS, X (KX, 1 ) ,  X (KX,  2 ) , 

C  1  TCEN, B, W, IER) 

C  THETA <KX)=Bf AVER 

C  EPS (KX) =W/VAR 

C  150  CONTINUE 
C 
c 
c 

SUBROUTINE  REMAV(PSI,M, AVER, SDV) 

C  ROUTINE  TO  CALCULATE  THF.  MEAN  AND  VARIANCE  OF  AN  ARRAY 

C  IT  ALSO  REMOVES  THE  MEAN  FROM  THE  ARRAY 

DIMENSION  PSIU) 

AVER=0 . 

SDV=0 . 

DO  10  1=1, M 

AVER= AVER  +  PS I ( I ) 

SDV=SDV+PSI (I) #*2 
10  CONTINUE 

AVER= AVER / FLOAT (M ) 

SDV=SDV/ FLOAT (M) -AVER  *  *  2 

IF  (M  .NE.  1)  SDV= ( FLOAT (M) /FLOAT (M-l ) )*SDV 
DO  20  1=1, M 

PSI ( I) =PSI ( I ) -AVER 
20  CONTINUE 
RETURN 
END 
C 
C 
c 

SUBROUTINE  SELECT) LIMIT, X, Y.TCEN, XOBS,T, PSI, 

1  TARSI , PARX, TOBS , H, MOPS , PI ST, TIM, a  I  fa , beta ) 

C  ROUTINE  TO  SELECT  THE  AT  MOST  “LIMIT*  NEARBY  POINTS 
C  TO  AN  INTERPOLATION  POINT  X,Y,TCEN 
C  LIMIT  IS  THE  MAXIMUM  NUMBER  OF  POINTS  TO  USE 
C  DIST  IS  THE  SPATIAL  RADIUS  OF  INFLUENTIAL  POINTS  IN  KM 

C  TIM  IS  THE  TEMPORAL  RADIUS  OF  INFLUENTIAL  POINTS  IN  DAYS 

DIMENS ION  XOBS ( 2 , 1 ) , T ( 1 ) , PS I { 1 ) . PARS 1(1), TOBS ( 1 ) 

DIMENSION  PARX (2,1) 

DIMENSION  INDEX { 2000 ), COR (2000) 
real  a,b 

COMMON / CB LOCK / C ( 20 ) 

EXTERNAL  F 
DATA  CPHSE/0 , 0/ 

NOBS=0 
DO  50  1*1, N 

DELX=X-XOBS (1,1) 

DELY= Y - XOBS ( 2 , I) 

DELT=TCEN-T ( I ) 

R=SQRT  ( ( DELX-CPHSE*  DEI.T)  *  *  2  ♦  DELY*  *  2 ) 

IF  (ABS(DELT)  .GT.  TIM)  GOTO  50 
IF  (R  .GT.  DIST)  GOTO  50 
NOBS = NOBS  «•  1 
INDEX (NOBS) = I 

COR { NOBS ) =F< DELX , DELY , DELT, *1  fa,  beta ) 


U  V 


50  CONTINUE 

IF  (NOBS  .EQ.  0)  GOTO  75 

IF  (NOBS  .GT.  LIMIT)  CALL  SORT (COR, INDEX, NOBS) 

IF  (NOBS  .GT.  LIMIT)  NOBS=LIMIT 
DO  70  1=1, NOBS 
>J=  INDEX  ( I ) 

FARX  < 1 , I ) =XOBS ( 1 ,  J ) 

PARX (2, I ) =  XOBS ( 2 , J ) 

TOBS(I)=T(J) 

FARSI ( I ) =PSI ( J ) 

C ( I ) =COR ( I ) 

70  CONTINUE 
75  CONTINUE 
RETURN 
END 
C 
C 
C 

SUBROUTINE  SORT (COR. INDEX , N) 

C  A  SHELL  SORT  ROUTINE  TO  SORT  INDEX  AND  COR  DOWN 

C  ACCORDING  TO  THE  VALUES  OF  COR 

DIMENSION  COR ( 1 ) , INDEX ( 1 ) 

IGAP=N 

5  IF  (IGAP  .LE.  1)  RETURN 
IGAP=IGAP/2 
IMAX=N- IGAP 
10  IEX=0 

DO  20  1=1 , IMAX 
I PLUSG= 1 4  IGA  P 

IF  ( COR ( I )  .GE.  COR ( IFLUSG) )  GOTO  20 
SAVE=COR ( I ) 

COR ( I ) =COR ( IFLUSG) 

COR ( I PLUSG ) =SAVE 
ISAVE= INDEX ( I ) 

INDEX ( I ) = INDEX ( I PLUSG ) 

INDEX ( I PLUSG) = ISAVE 
IEX=1 
20  CONTINUE 

IF  (IEX  .NE.  0)  GOTO  10 
GOTO  5 
END 
C 
C 
C 

SUBROUTINE  OBJAN( PSI , L, T, N, X, Y, TCEN, B,W, IER, alfa, beta) 
C  THE  SPACE-TIME  OBJECTIVE  ANALYSIS  ROUTINE 
C  VERSION  FOR  1  INTERPOLATION  FOINT 
USES  2  SPACE  AND  1  TIME  DIMENSION 
NOTE  DELTA  T  =  TCEN  -  T(J) 

C  L  IS  THE  ARRAY  OF  OBSERVATION  POSITIONS,  IN  KM 

C  T  IS  THE  TIME  OF  OBSERVATION,  IN  DAYS 

C  X  IS  THE  ARRAY  OF  INTERPOLATION  POSITIONS,  IN  KM 

C  TCEN  IS  THE  CENTRAL  INTERPOLATION  TIME 

C  PSI  IS  THE  ARRAY  OF  OBSERVATION  VALUES 


n  n 


C  B  IS  THE  INTERPOLATED  VALUE 

C  W  IS  THE  INTERPOLATION  ERROR  LIMIT 

C  N  IS  THE  NUMBER  OF  OBSERVATION  POINTS 

C  IER  IS  AN  ERROR  FLAG,  ZERO  FOR  NO  ERROR 

C  -1  No  data  (WARNING) 

C  -3  Poor  matrix  inversion  (WARNING) 

DIMENSION  PSI(l) ,T(1) ,L(2, 1) 

COMMON/CBLOCK /C ( 2  0 ) 

REAL* 8  A(20, 20) 

REAL  L 

COMMON/ ERR /E 
EXTERNAL  F 
B=0 . 

W=1.0 
IER=- 1 

IF  (N  .LE.  0)  GOTO  500 

C  CALCULATE  THE  INVERTED  AUTOCORRELATION  MATRIX  FOR  THE  OBSERVATIONS 

CALL  SETA(A,L,T,N, IER, al fa, beta) 

IF  (IER  -GT.  0)  GOTO  500 
CALL  GETAVE(A, PSI , M, AVE) 

C  CALCULATE  THE  MATRIX  C 

C  --  already  calculated  in  this  version,  common  block  CBLOCK 
C 

C  CALCULATE  THE  RMS  INTERPOLATION  ERROR, W 
C  AND  CALCULATE  THE  INTERPOLATED  VALUE  B 
W=0. 

B=0 . 0 

DO  150  1*1, N 
H=0 . 0 
DUMC=C ( I ) 

DO  140  J  = 1 , N 

P=DUMC*C ( J ) *  SNGL ( A ( I , J ) ) 

W=W»P 

P2=SNGL (A ( I , J ) ) *  PSI ( J) 

H=IHP2 

140  CONTINUE 

DUMY*DUMC*H 
B=B+DUMY 
150  CONTINUE 
B=B»  AVE 
W^ABS(l.-W) 

500  CONTINUE 
RETURN 
END 
C 

c 

c 

SUBROUTINE  SETA (A, PARX, T, NOBS, IER, alfa, beta) 

C  THIS  ROUTINE  CALCULATES  THE  AUTOCORRELATION  MATRIX  FOR  THE 
OBSERVATIONS  GIVEN  THE  POSITIONS,  PARX  AND  TIMES,  T 
IT  RETURNS  THE  INVERTED  MATRIX 
DIMENSION  PARX ( 2 , 1 ) . T ( 1 ) 

REAL* 8  A(20, 20) , Det 
Integer  IP(40) 


COMMON/ERR/E 
EXTERNAL  F 
DATA  NA/20/ 

C  The  Guard  value  for  DETERMINANT  WARNINGS 
DATA  GUARD/ 1 . OE-4/ 

1  FORMAT  (5X, 'MATRIX  A  IS  SINGULAR') 

2  FORMAT  (5X, 'ERROR, MATRIX  A  IS  TOO  SMALL',/, 

X  '  NA  MUST  BE  .GE.  NOBS',/,'  NA= ' , 13 , 5X, 'NOBS= ' , 13 , // ) 

3  FORMAT  (5X, 'WARNING,  DETERMINANT  IS  VERY  SMALL  ( * , 1PE11. 4, ' ) ' , 

X  ’  —  TRY  SMALLER  NUMBER  OF  INFLUENTIAL  POINTS') 

C  TEST  THE  SIZE  OF  THE  OBSERVATION  ARRAY 
IER=1 

IF  (HA  .LT.  NOBS)  PRINT  2,NA,NOBS 

IF  (NA  .LT.  HOBS)  RETURN 

IER=0 

DO  20  1=1, NOBS 
DO  10  J = I , NOBS 
DELT=T( I ) -T( J) 

DELX=PARX  (1,1)-  PARX  ( 1 ,  •! ) 

DELY=PARX  ( 2 ,  I ) -PARX  ( 2 ,  .7 ) 

A ( I , J) =DBLE(F (DELX, DELY, DELT, alf a, beta) ) 

A(J, I ) =A ( I . J ) 

10  CONTINUE 

A ( I , I ) = A ( I , I) +  DBLE ( E ) 

20  CONTINUE 

C  INVERT  THE  NOBS* NOBS  MATRIX  A 

C  THE  INVERTED  MATRIX  IS  NAMED  A 

Call  I nvMt x { A , NA , A , NA , HOBS , De t , IP, Ier, alfa, beta) 

IF  (IER  .HE.  0)  PRINT  1 
IF  (IER  .NE.  0)  GOTO  40 
C  CHECK  THE  DETERMINANT 

IF  (DET  .LT.  GUARD)  PRINT  3 , DET 
IF  (DET  .LT.  GUARD)  IER=-3 
40  CONTINUE 
RETURN 
END 
C 
C 

c 

SUBROUTINE  GETAVE(A,  PS  I,  II,  AVE) 

C  Calculate  and  remove  the  weighted  mean 

DIMENSION  FSI ( 1 ) 

DIMENSION  C(20),D(20) 

REAL* 8  A(20, 20) 

DO  10  1=1, N 
C(I)=0 
D(  I ) =0 
DO  10  K=1,N 

C(I)=C(I)4A(I,K)*PSI(K) 

D(I)=D(I)4A(I,K) 

10  ENDDO 
SUM1=0 
StJM2=0 
DO  20  1=1, N 


SUM1=SUM1 +C ( I ) 
SUM2=SUM2+D( I ) 

20  ENDDO 

AVE=SUM1/SUM2 
DO  30  1  =  1, N 

PS  I  ( I )  =  PS  I (I)- A VE 

30  ENDDO 
RETURN 
END 


c  APPENDIX  6B 


C  tliia  is  the  program  that  links  all  th<?  elements  together. 
r  and  also  controls  the  reduction  or  removal  of  successive  bathys 

c 

parameter (most  =  56) 

tea  1  alfa(4),beta(4), amode (most ) , xpos (most ) , ypos (most ) 
real  a,b 

integer  n,pr inter , loop, ran (most ), i , j , left , ans, num  .reply 

c  read  in  the  data  for  each  time  data  set  is  reduced  by  one 
c 

do  100  2  =  l.most 

write!*,*)  ’  process  1  =  y  2  =  end' 
read (  *  ,  *  )  reply 
if (reply. eq. 2)  goto  200 

write!  *,*)  'how  many  XBTs  do  you  want  to  use  (  max  133) 
read ( * , * )  ans 

j  =  most  -  at’.s 

call  reduce! j) 
left  =  most  -  j 

c 

c  read  in  modes 


do  10 

n  =  1,4 

open 

(  unit  = 

1, 

file  = 

’  ldecpos .mat ' 

open 

(  tin  it  = 

11. 

file  - 

’ rmodl ' ) 

Open 

(  unit  = 

12, 

file  = 

' rmod2 '  ) 

open 

(  unit  = 

13, 

file  = 

' t  mod  3  ' ) 

open 

(  unit  = 

14, 

file  ^ 

'  rt»od4  ’ ) 

c 

c 

c  set  parameters  for  input 
c 

c  chose  values  for  a  and  b 
c 

alfa(l)  =  155.5 
beta ( 1 )  =  85.3 

alfa(2)  =  92.15 

beta (2)  =  62.14 

alfa(?)  =  36.65 

beta ( 3 )  =  28.35 

a  1  fa  ( 4 )  =  39.4 

beta (4)  =  30.55 

c 

c 

C 


c  now  loop  through  oa  four  times. once  for  each  mode 
c 


c  read  in  latitude  and  longitude  of  obs 

do  26  i  =  1 ,  1  e  f  t 

read(l,*)  num, xpos ( i ) , ypos ( i ) 

26  continue 
rewind ( 1 ) 
c 

c  read  in  modes 

do  27  1  .  1 , left 

read (10  +  n,*)  amode(i) 

2‘,  continue 

close ( 1 ) 
close (11) 
close ( 12 ) 
close ( 13 ) 
close( 14) 
c 

c  give  values  of  a  and  b 
c 

a  =  alfa (n) 
b  =  beta (n) 

c  set  output  diagnostics, 
c 

printer  =  30  +  n 
loop  =  n 

c 

c 

call  modeoa (amode, xpos , ypos, a , b, printer , lef t , loop) 


10  continue 

c  print  out  modes  and  errors  to  combined  files 

call  allmods 
100  continue 

200  write(*,*)  'program  finished’ 

end 


c 


Program  by  E.F.  CARTER 


subrout ine  modeoa (amode, x, y, al fa,  beta ,  printer , most , loop) 
c  PROGRAM  modeoa 

C 

c 

c 

C  UNIT  1  IS  THE  XBT  (INPUT)  DATA 

C  UNIT  2  IS  THE  PRINTABLE  OUTPUT  DATA  (DIAGNOSTICS) 

C  UNIT  4  IS  THE  UNFORMATTED  OUTPUT  DATA 

C 

c  INTEGER  INFILE, DISK, PRINTER 

c  PARAMETER  (INFILE=1,  DISK=4) 

C 

Parameter  (MXOBS  =  56) 

INTEGER  DAY (MXOBS) ,Gmt (MXOBS) , pr inter , most 
DIMENSION  X (MXOBS) ,Y (MXOBS) .Tin (MXOBS) , amocle (MXOBS) 
DIMENSION  XOBS ( 2 , MXOBS ) , TOBS (MXOBS) , Ur BS (MXOBS) 

DIMENSION  UOPT(20) , VOPT ( 20 ) , TOPT( 20 ) , XOPT(2, 20) 

DIMENSION  XI(119) , YI(119) ,UI(119) , ERRU{119) 

INTEGER  START, EMPTY, loop 
Real  MnL ay, MxDay, al fa, beta 
COMMON/ERR/E 
EXTERNAL  F 
DATA  EMPTY/ 156/ 
c  DATA  MOST/ 57/ 

DATA  LIMIT/ 5/ 

C  SPATIAL  AND  TEMPORAL  LIMITS 
DATA  DIST, TIM/ 150 . , 0 . / 

DATA  CFX.CFY/ 15.0, 15.0/ 

DATA  START/ 1/ 

DATA  TINC/1/ 

C 

DATA  NOBJ/1/ 

r* 

£  **4++*4*+4****++****4****4****#****4+«*+4+ft**+****+4*«**«* 

c  write (*,*)  'modeoa  now  being  called' 

c  open  files  Cor  output 


open 

(unit 

= 

31. 

file 

= 

' out put  1' ) 

open 

(unit 

= 

41. 

file 

- 

'gridmod . 1 ’ ) 

open 

(unit 

= 

51. 

file 

= 

'errormod. 1 ' ) 

open 

(unit 

= 

32, 

file 

'output  2 ' ) 

open 

(unit 

= 

42, 

file 

- 

'gridmod. 2 ’ ) 

open 

( un  i  t 

= 

52, 

file 

- 

'errormod. 2 ' ) 

open 

(unit 

= 

33, 

file 

'output  3 ' ) 

open 

(unit 

= 

43. 

file 

= 

'gridmod. 3 ’ ) 

open 

( un  i  t 

= 

51, 

file 

a 

'errormod. 3 ' ) 

open 

(unit 

= 

34, 

file 

- 

'output  4 ' ) 

open 

(unit 

= 

44. 

file 

tr 

' gr idmod . 4 ' ) 

open 

(unit 

= 

54, 

file 

= 

’ errormod. 4 ' ) 

c 


c 


1  FORMAT 

2  Format 

3  Format 

22  Format 

23  Format 

4  FORMAT 

5  FORMAT 

6  FORMAT 
& 

7  FORMAT 

8  FORMAT 

9  FORMAT 
10  Format 
18  FORMAT 


(5X, 'DIAGNOSTICS  OF  ALL  modal  OBSERVATIONS 
(5x, 'X  Position  Diagnostics  (Km):  ') 

(5x, 'Y  Position  Diagnostics  (Km):  ’) 

(5x, 'X  grid  Diagnostics  (Km):  ') 

(5x,  *Y  grid  Diagnostics  (Km):  ') 

(5X, ' INTERPOLATED  mode  FIELD  DIAGNOSTICS:’) 

(5X, ' INTERPOLATED  mode  Error  FIELD  DIAGNOSTICS:') 
(5X, 'JULIAN  DATE: '.F8.3, 

'  NUMBER  OF  Observations  used  :',I4) 

(5X,  'Minimum  date:'F8.3,'  maximum  date:\f8.3) 
(7X, 'ERROR,  XOBS  TOO  SMALL.  MXOBS=',I3) 

(5X, 'Assumed  NOISE  LEVEL: ' , F8 . 5) 

(5X, 'Date  (Time)  Diagnostics  (Julian  Date):  ') 
(5X, 'Humber  of  influential  points  is  :  ',14,//) 


£,  ***************************************************  ********** 


E=0 . 1 

WRITE  (PRINTER, 9)  E 

c 

£.  ************************************************************** 


C  READ  IN  THE  OBSERVATION  DATA 
c  1  =  1 

c  25  Continue 

c  OPEN  (UMIT= INFILE, Fi le= ' xbt . dat ' ) 

c  READ  ( INFILE, * , End=30)  Day ( I ) ,Cmt ( I ) , Y ( I ) , X ( I ) , SST( I ) 

c  1=1+1 

c  Goto  25 

c  30  Continue 

c  Most=I-l 

c  CLOSE  (UHIT= INFILE) 

r  *************************************************************** 

c 

c  Read  in  the  observed  data. 

c  for  day  and  time 

do  25  i  =  l,most 
Day  ( i )  =  1 

Gmt(i)  =  1 

tin(i)  =  1 

25  continue 

c  for  latitude  and  longitude  of  observations 
c 

c  open  (unit  =  I,  file  =  'dec.pos') 

do  26  i  =  l,most 

write!*,*)  x(i),y(i) 

26  continue 
c 

c  for  modal  amplitudes 
c 

c  open  (  unit  =  2,  file  =  'xaa') 


u  u 


c  do  27  i  =  l.most 

c  read (2,*)  amode(i) 

c 27  continue 

c 

c  for  latitude  and  longitude  of  grid  points 
c 

open(unit  *  3,  file  a  'grid.poa') 
do  28  i  =  1,119 

read(3, * )  n , xi ( i ) , y i ( i ) 

28  continue 

rewind(3) 
close ( 3 ) 

C  **************♦*******<***  **********************+**<l**il******A 


CALL  SCALE (X, Y,Most) 


Write  (Printer, 2) 

Call  Diag(X, Most , Printer ) 
Write  (Printer, 3 ) 

Call  Diag(Y, Most, Printer) 
c 

CALL  SCALE(Xi,Yi, 119) 
c 

Write  (Printer, 22) 

Call  Diag(Xi, Most, Printer) 
Write  (Printer, 23' 

Call  Diag(Yi, Most, Printer) 

c 

c  CALL  JULIAN ( DAY, Gmt, Tin, Most) 
c 

Write  (Printer, 10) 

Call  Diag(Tin, Most, Printer ) 


WRITE  (PRINTER, 1) 

CALL  DIAG (amode, MOST, PRINTER ) 
c 

Call  Remav(amode, Most , Ave, V) 
c 

r  +***+««*4****4**4****«4**+4****4*4««*«*****ft*«*4*4**#*4 


c 

c 

c 

c 

c 

c 

c 

c 

r 

c 


SET  UP  THE  INTERPOLATION  POSITIONS 
M-0 

DO  40  * J  —  1 #  21 
DO  40  1  =  1,  33 
M=M*1 

XI (M) =CFX* ( I- 17 ) 

YI (M) =CFY# (J-ll) 

40  ENDDO 


DO  SEVERAL  ANALYSES 


WRITE  (PRINTER, 18)  LIMIT 
TCEN=START-TINC 
DO  500  IOBJ=l , NOBJ 
TCEN=TCEN+TINC 

C  GET  THE  USABLE  OBSERVATIONS  FOR  THIS  DATE 

CALL  GETOBS (Tin, X, Y,amode, Most, TCen, Tim, 

X  Tobs , Xobs , UOBS ,  N) 

MNDAY=TCen-TIM 
MXDAY=TCen  *-TIM 
WRITE  (PRINTER, 6)  Tcen.N 
WRITE  (PRINTER, 7)  MNDAY, MXDAY 
WRITE  (PRINTER, 18}  LIMIT 
IF  (N  .EQ.  0)  GOTO  500 
IF  (N  .GT.  MXOBS)  THEM 

WRITE  (PRINTER, 9)  MXOBS 
GOTO  500 

EHDIF 

C 

C  m  is  number  o£  grid  positions, and  will  not  change, 
m  =  119 
DO  50  k  si,  m 

CALL  SELECT ( LIMIT, X I ( K ) , Y I ( K ) , TCEM , XOBS . TOBS , UOBS , 
X  UOFT, XorT, TOPT, N, NOBS, DIST, TIM, al £a, beta ) 

C  Print  *, 'M,  Nobs,  ',N,N obs 

C  Do  49  I PXJ  = 1 , NOBS 

C  Print  * , Xopt (1,1 PXJ ) , XOpt (2,1 PXJ ) . Uopt ( I PXJ ) 

C  49  EndDo 

CALL  OB.!AN(UOPT.  XOPT,  TOPT,  NOBS,  XI  (K)  ,YI(K)  , 

X  TCEN, UI  (K) , ERRU(K) , IER, alfa, beta) 

UI (K) =UI (K) +Ave 
50  ENDDO 

c 

WRITE  (PRINTER, 4 ) 

CALL  DIAG(UI,M, PRINTER) 

WRITE  (PRINTER, 5) 

CALL  DIAG(ERRU,M, PRINTER) 


WRITE 

(30 

4 

loop, * 

)  'centre  time’.TCen 

WRITE 

(30 

4 

loop, * 

)  '  number 

of  obs  point 

WRITE 

(50 

4 

loop, * 

)  -72.  -64, 

37  .40 

write 

(50 

4 

loop, * 

)  17,  7 

WRITE 

(40 

4 

loop, * 

)  -72,  -64, 

37  ,40 

write 

(40 

4 

loop  , 

♦)  17.  7 

• 

do  51 

i  = 

1. 

m 

WRITE 

(40 

4  loop, * ) 

UX  (it 

WRITE 

(50 

4  loop,  * > 

ERRU(i) 

51  continue 


n  o 


500  CONTINUE 

close (31) 
close ( 32 ) 
close<33 ) 
close ( 34 ) 
close(41) 
close ( 42 ) 
close (43) 
close (44 ) 
close (51 ) 
clone ( 52 ) 
close(53 ) 
close(54) 

END 

£  ******************************************************** 

C  ******************************************************** 
f-.  ********************************************************* 


SUBROUTINE  GETOBS ( Day , Posx , Posy , UOBS, Ninp, CDay , Tim, T. X, U, N) 
C 

C  Input  data : 

C  Day , Posx, Posy 

C  UOBS 

C  Ml  tip 

C  CDay 

C  TIM 

C 

C  Output  data: 

C  T,X  location  of  data 

C  U  chosen  observation  data 

C  N  number  of  points  used 

C 

ROUTINE  TO  GET  THE  OBSERVATION  DATA 
N  IS  THE  NUMBER  OF  OBSERVATIONS 
C 

DIMENSION  T(1),X(2,1),U(1) 

DIMENSION  Day ( 1 ) . UOBS ( 1 ) . POSX ( 1 ) , POSY ( 1 ) 

Real  MXDAY , MNDAY 
MXDAY=CDAY4TIM 
MNDAY- CDAY-TIM 
KOUNT* 0 

Do  10  1=1, Ninp 

IF  (DAY (I)  .GT.  MXDAY)  GOTO  10 
IF  ;DAY(I)  .LT.  MNDAY)  GOTO  10 
KOUNT=KOWmi 
T(KOUNT) =DAY ( I ) 

X ( 1 , KOUNT) =POSX ( I ) 

X ( 2 , KOUNT) =POSY ( I ) 

U (KOUNT) =UOBS ( I ) 

10  Continue 
N- KOUNT 
RETURN 


space-time  location  of  data 
observed  data 
number  of  input  points 
Central  day  of  estimate 
width  of  time  window 


on  o  o  o 


END 

C 


SUBROUTINE  SCALE(X,Y,N) 

Scale  Lat  and  Long  to  Km 
DIMENSION  X(l)  .  Y<1) 

Parameter  (XCen  =  -70.  S,  YCen  -  3R.25) 
Parameter  (CFX  *  110.99,  CFY  =  87.84) 
DO  10  1=1, N 

Y(I)=CFYM  Y(I)-YCen) 

X(I)=CFX*(  X(I)-XCEN) 

10  EMDDO 
RETURN 
END 


C 

SUBROUTINE  JULIAN ( DAY , Gmt , Time , N ) 

INTEGER  DAY ( 1 ) , Gmt ( 1 ) 

Dimension  Timed) 

Integer  Offset 

Parameter  (Offset  =  86000,  JulianBS  =  6431) 

Real  MnDay.MxDny, Minutes 

1  FORMAT  (5X. 'MIN  DATE  : ' , F8 . 3 , '  MAX  DATE  : ' , F8 . 3 
MNOAY=9999 . 9 
MXDAY =0.0 
DO  10  1=1, N 

Time ( I ) =Float (Gmt ( I ) ) / 100 . 

Minutes=  100 .  *  (Time  ( I )  -  IFix(Timed) ) ) 

Time ( I) =Time ( I )  +  Minutes/60.0 

Timed)  =  Float  (Jul  ianR6  f  Day(I)  -  Offset) 

IF  (Time (I)  .GT.  MXDAY)  MXDAY =Time ( I ) 

IF  (Time(I)  .LT.  MNDAY)  MNDAY=Time ( I ) 

10  ENDDO 

WRITE  (2,1)  MNDAY, MXDAY 
RETURN 
END 
C 
C 

c 

FUNCTION  F ( X, Y , T, a 1 f a , beta ) 

C  THE  CORRELATION  FUNCTION 

C  THE  SCALE  FACTORS 

c  Parameter  (a=111.6,  b=B6.6) 

a  =  alfa 

b  =  beta 

r2=x*  *2  +  y**2 

F= (1 . 0  -  r2/a**2) *exp(-r2/b**2) 

RETURN 

END 


Time (I ) /24 . 0 


C  this  program  generates  first  158  numbers  ranomly 
c 

parameter (max  =  58) 
integer  i,m,x(max),n 

open (  unit  =1,  file  =  'randnumdeep' ) 
call  srand (3) 
do  10  i  =  1  ,  max 
*0  n  =  irandO 

x(i)  =  n 

if(i.gt.l)  then 

do  20  m  =  1 , i- 1 

if  ( (x ( i ) . eg . x (m) ) .or . (x ( i ) . gt .max-2) )  goto  40 
20  continue 

end  if 


write ( * , 1 )  i , x ( i ) 
write(l,*)  x(i) 

10  continue 
end 


close  < 1 ) 
close (2) 
close (11) 
close(12) 
close( 13 ) 
close(14) 
close{20) 
close(21) 
close (22) 
close(23 ) 
close(24) 
end 


o  o 


c  APPENDIX  7B 

C  Program  by  E.F.  Carter 

SUBROUTINE  INVMTX  (A, NA , V, NV, N, D, IP, IER.  al £a ,  beta ) 

C  Double  Precision  version 
C  MATRIX  INVERSION  V=INV(A) 

C  THE  ARRAY  A  MAY  BE  ENTERED  AS  V  TO  SAVE  MEMORY 
C  IP  MUST  BE  DIMENSIONED  TO  AT  LEAST  2‘N 
INTEGER  NA, NV. N. IP(t ) , !ER 

REAL* 8  A (HA, N)  ,V(NV,N)  ,D 

Real *8  VMax,  VH,  PVT,  PVTMX,  HOLD 

I EXMAX  IS  SET  TO  THE  LARGEST  BASE  TEN  EXPONENT  THAT  CAN  BE 
REPRESENTED  ON  THE  MACHINE,  I.E.  LARGEST=10*  *  I EXMAX 
DATA  I EXMAX/38/ 

115  FORMAT ( 28H0 *MATR IX  SINGULAR  IN  INVMTX') 

116  FORMAT ( 3 4H0* DETERMINANT  TOO  LARGE  IN  INVMTX*) 
tER  =  IERINV(N, NA, NV) 

IF  (IER  -HE.  0)  RETURN 
DO  102  J=1 , N 

IP(J)  =  0 
DO  101  1=1, M 

V(I,J)  =  MI.J) 

101  CONTINUE 

102  CONTINUE 
D  =  1. 

I  EX  =  0 
DO  110  M=1,N 

VMAX  =  0. 

DO  104  J*1,M 

IF  (IP(J)  .ME.  0)  GO  TO  104 
DO  103  1=1, N 

IF  (IF(I)  .HE.  0)  GO  TO  103 
VH  =  ABS (V ( I , J ) ) 

IF  (VMAX  .GE.  VH)  GO  TO  103 
VMAX  =  VH 
K  =  I 
L  =  J 

103  CONTINUE 

104  CONTINUE 
IP(L)  =  K 
NFM  =  N+M 
IP(NPM)  =  L 
D  =  D*V (K, L) 

105  IF  (ABS (D)  .LE.  1.0)  GO  TO  106 

D  =  D*0. 1 
I EX  =  IEX+1 
GO  TO  105 

106  CONTINUE 
PVT  =  V(K,  L) 

IF  (M  .EQ.  1)  PVTMX  =  ABS (PVT) 

IF  ( ABS < PVT/FLOAT (M) )+ PVTMX  .EQ.  PV1MX)  GOTO  113 
V(K,L)  =  1. 

DO  107  J=1,N 


subroutine  reduced ) 

parameter (most  =  56) 

integer  i , j , ran (most ), va  1 , p 

real  amodl (most ) , amod2 (most ) , a mod 3 (most ) 

real  amod4 (most ) , xpos (most) . ypos (most ) 

c  read  in  original  modes 


open  (unit  =  2,  file  =  ' randnumdeep' ) 


open 

( 

unit 

= 

1, 

file 

open 

( 

unit 

= 

11. 

file 

open 

( 

unit 

= 

12, 

file 

open 

( 

unit 

= 

13. 

file 

open 

( 

unit 

= 

14. 

file 

c 

c  write  out 

reduced 

set . 

open 

( 

unit 

= 

20, 

file 

open 

( 

unit 

- 

2 1  > 

file 

open 

( 

unit 

= 

22, 

file 

open 

( 

un  i  t 

= 

23, 

fils 

open 

( 

unit 

= 

24, 

file 

c 

r.  **************************** 


’  decpos . mat ' ) 
'  mod 1 . ma  t ' ) 
’mod2 .  mat. '  ) 
'modi . mat ' ) 
'mod4 .mat ' ) 


' rdecpos .mat ' ) 

' rmodl ’ ) 

’ rmod2 ’ ) 

' rmodl ' ) 

' rmod4 '  ) 

**************************** 

by  j,the  parameter  fed  from 


c 

c  this  section  will  reduce  data  set 
c  program  driver. 


p  =0 

do  100  i  =  l,most 

read (2, *)  ran  I  i) 
write (*,*)  i,  tan(i) 
read(l,*)  xpos(i) ,ypos(i) 
road (11, * )  amodl ( i ) 
read (12,*)  amod2(i) 
read (13,*)  amod3(i) 
read (14,*)  amod4(i) 

100  continue 


do  110  i  =  l.most 
do  120  n  =  1 , j 

if ( i . ne . ran (n) )  goto  120 
i f ( i . eg. ran (n ) )  goto  110 
120  continue 

P  =  P  ♦  1 

write(20,*)  p,  xpos  ( i ) ,  ypos  (  i)_ 
write(21,*)  amodl(i) 
write(22,*)  amod2 ( i ) 
write(23,*)  amod3(i) 
write (24,*)  amod4( i ) 

110  continue 


MOLD  =  V(K,J> 

V(K,J)  =  V(L,J) 

V  ( L,  J )  =  HOLD/ PVT 

107  CONTINUE 

DO  109  1=1, N 

IF  (I  .EQ.  L)  GO  TO  109 
HOLD  =  V(I,L) 

V( I , L)  =  0. 

DO  108  J=1 , N 

V(I,.J)  =  V(I, J)-V(L, J) ‘HOLD 

108  CONTINUE 

109  CONTINUE 

110  CONTINUE 
M  =  N*N+ 1 

DO  112  J=1 , N 
M  =  M-l 
L  =  IP (M) 

K  =  IP(L) 

IF  (K  .EQ.  L)  GO  TO  112 
D  =  -D 

DO  111  1=1, N 

MOLD  =  V( I , L) 

V( 1 , L)  =  V(I,K) 

V(I,K>  =  HOLD 

111  CONTINUE 

112  CONTINUE 

IF  ( IEX  .GT.  I EXMAX)  GO  TO  114 

D  =  D* 10 .  “  IEX 

RETURN 

113  IER  =  33 
PRINT  115 
RETURN 

114  IER  =  1 

D  =  FLOAT (IEX) 

PRINT  116 

RETURN 

END 

FUNCTION  IERINV  (N.HA.MV) 

103  FORMAT (23H0*  N  .LT.  1  IN  IMVMTX  *) 

104  FORMAT(24MO*  NA  .LT.  N  IN  INVMTX  *) 

105  FORMAT (24 M0*  MV  .LT.  M  IN  INVMTX  *) 
IERINV  =  0 

IF  (M  .GE.  1)  GO  TO  101 
IERINV  =34 
PRINT  103 
RETURN 

101  IF  (NA  .GE.  N)  GO  TO  102 
IERINV  =  35 

PRINT  104 
RETURN 

102  IF  (NV  .GE.  N)  RETURN 
IERINV  =36 

PRINT  105 
RETURN 


%  APPENDIX  8B 


'  tNig  matlab  file  reconstruci.r,  the  bathys  at  the  qrid  points 
'  using  4  modes  only,  after  the  OA  with  reducing  number  of 
'  initial  XBT's. 

% 

%  this  part  reconstructs  all  synthetic,  grid,  XBTs  using  4  modes. 

' 

clear 

clg 

hold  off 


load  vec; 
load  gridmod; 
load  errmod; 
grid  =  gridmod'; 

' 

for  j  =  1:4; 
for  i  =  1:80; 

eigvec(i,j)  =  vec ( i ,  ( 8 1- j  )  )  ; 

end 

end 


recbath  -  eigvec*gr id; 

save  recbath. mat  recbath  /ascii; 

\  PLOT  GIVEN  XBT 

'  this  part  plots  a  given  ha  thy,  .SELECT  XBT  NUMBER 
'  1,  being  lower  left  corner. 


'  maxis  =  JO  30  -80  0 } ; 

'  axis (maxis); 

'plot (recbath(56) ,'.'); 

'  grid; 

't it le ( '  XBT  56  '  )  ; 

'ylabel  ('depth  metres'); 

'  xlabeK'temp  Celcius'); 

'  hold  on 
'  pause 
'print 

i  ummuummmnmmmmmmmmmmunm 

' 

1  now  to  calculate  the  error  in  each  reconstructed  bathy 
eigvecsrd  a  eigvec’  *  eigvec; 
xbtvar  =  eigvecsrd*errmod* ; 
xbtvar  =  sum (xbtvar) 
view  =  reshape (xbtvar, 17,7) 


view  =  view* • 


view  =  flipud  (view); 

<view  s  nerr; 

save  batherr.tnat  view  /ascii; 


hold  off 

%  plot  routine  for  tlATLAB,  of  the  reconstructed  oa  error. 

for  i  =  1:17; 

x(i)  =  -72.5  +0 . 5* i ; 

end 

for  i  =  1:7; 

y(i)  =  36.5  +  0 . 5  * i ; 

end 

c  a  contour (view, x, y ) ; 
v  =  (0. 2, 0.4, 0.6, 0.8); 
clahel (c, v) ; 
grid; 

%  t it le (' reconstruct  ion  error,  using  1  out  of  156  XBTs');  ^  the  title  must  chn 
xlabel ( 'Longitude  west  from  Greenwich');  *  each  time 

ylabel ( 'Latitude  North'); 

title('  reconsruction  error  variance  using  all  156  XBTs') 
hold  on 


%  mmtmmmmmmmimmmmmmmn 

% 

% 

%  loading  in  data  positions. 

load  rdecpos 

long  =  rdecpos) :, 2) ; 
lat  =  rdecpos) : , 3 ) ; 

^plot ( long, lat 

Ixlabel (' Longitude  west  from  Greenwich') 
lylabel { ’Latitude  North') 

^title  ('position  of  deep  XBTs') 

^text (-70.5,38.5, 'X' ) 

Uprint 

tpause 

nummmnmmmmnmmmmt 

%  mark  on  box 
box  =(  -70.5  38.1 


-70.5  38.75 

-69.25  38.75 

-69.25  38.5 

-70.5  38.11 

blong  =  box( : , 1 ) 
blat  =  box( : ,2) 
tplot (blong, blat ) 
pause 
print 

mmtummmmmumtmmtmmt 

%  load  in  each  modal  error. 

hold  off 
clg 

moderrl  a  errmod ( * , 1 ) ' ; 

moderrl  =  reshape (moderr 1 , 17 , 7 ) ; 

moderrl  =  moderrl'; 

moderrl  =  f 1 ipud (moderr 1 ) ; 

kl  =  contour (moaerr 1 , x, y ) ; 

clabel (kl, v) ; 

title  ('mode  1  error') 

Sprint 

%pause 

moderr2  =  errmod {: ,2) ' ; 

moderr2  =  reshape (moderr2 , 17 , 7 ) ; 

moderr2  =  moderr2'; 

moderr2  =  f 1 ipud (moderr2 ) ; 

k2  =  contour (moderr2, x, y) ; 

clabel (k2, v) ; 

title  (’mode  2  error') 

Iprint 

Ipause 

moderr3  =  errmod (:, 3 )' . 

moderr3  *  reshape (moderr 3 , 17 , 7) ; 

moderr3  =  moderr3'; 

moderrB  =  Elipud(moderr3 ) ; 

k3  =  contour (moderr 3 , x, y) ; 

clabel (k3,v) ; 

title  ('mode  3  error') 

%pr int 
fpause 

moderrl  =  errmod ( : , 1 ) ' ; 

moderr4  =  reshape (moderr 4 , 17 , 7 ) ; 

moderrl  =  moderr4 ' ; 

moderrl  =  f 1 ipud (moderr 4 ) ; 

kl  =  contour (moderr4,x,y) ; 

clabel (kl, v) ; 

title  ('mode  4  error') 

tprint 

%pauae 


subplot (221) , contour (moderr 1, x. y )  ; 
subplot  (222) ,  contour  (moderr2,x,y)  ; 
subplot (223) , contour (moderr 3, x,y) ; 
subplot (224) , contour (moderr4, x,y) ; 
Sprint 
Dclg 


ummmmmimimmmmmmn 

%  plot  routine  for  CONTOUR 
Dnerr=  nerr'; 

tsave  err. con  nerr  /ascii; 


t 

* 

tpluserr  =  recbath( : , 21)  +  totalerr ( : . 21) ; 
Dminuserr  =  recbath( : , 21 )  -  totalerr { : . 21) ; 
Dplot  (pluserr) 

Dplot (minuserr) 

Dglunk  =  totalerr (:, 21) ; 

Delong  =  reebath ( t , 21 ) ; 

Dsave  err7105.mat  glunk  /ascii 
isave  bath7105.mat  clong  /ascii 
ifor  j  =  1  :  1 
Ifor  i  =  1:80 

Ddeptherr (i. j)  =  glunk ( i , j )* 100/clong ( i . j ) ; 

D»nd 

lend 

Dsave  percenerr .mat  deptherr  /ascii 
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